{
"cells": [
{
"cell_type": "markdown",
"id": "f53f2418",
"metadata": {},
"source": [
"# Open Method\n",
"**강좌**: *수치해석*"
]
},
{
"cell_type": "markdown",
"id": "24713c71",
"metadata": {},
"source": [
"## Newton-Raphson Method\n",
"\n",
"### 알고리즘\n",
"* 함수 $f(x)$ 가 미분 가능하고 연속함수인 경우에 대해서 Taylor Expansion 을 이용하면 다음과 같다.\n",
"\n",
"$$\n",
"f(x_{i+1}) = f(x_{i}) + f'(x_i)(x_{i+1} - x_i) + \\frac{f''(\\xi)}{2!} (x_{i+1} - x_i)^2.\n",
"$$\n",
"\n",
"여기서 $\\xi \\in (x_i, x_{i+1})$ 이다. 선형 근사하고, $x_{i+1}$ 에서 x 축과 교점이 생기는 경우 다음과 같다.\n",
"\n",
"$$\n",
"0 = f(x_{i+1}) = f(x_{i}) + f'(x_i)(x_{i+1} - x_i).\n",
"$$\n",
"\n",
"즉, 아래 식을 반복적으로 해석해서 계산할 수 있다.\n",
"\n",
"$$\n",
"x_{i+1} = x_{i} - \\frac{f(x_i)}{f'(x_i)}\n",
"$$\n",
"\n",
":::{figure-md} newton-fig\n",
"
\n",
"\n",
"Newthon Method (from Wikipedia)\n",
":::\n",
"\n",
"#### 종료 판정 기준\n",
"Bracketing method와 동일하게 근사 상대 오차 $\\epsilon_a$ 의 크기가 $tol$ 보다 작으면 멈춘다.\n",
" \n",
":::{note}\n",
"$\\epsilon_a< tol$ 이면 멈춘다.\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "7a9f3812",
"metadata": {},
"source": [
"### 예제 1\n",
"$f(x) = e^{-x} - x$ 의 근을 구하시오. 초기 가정은 $x_0 = 0$ 이다.\n",
"\n",
"#### By Hand\n",
"도함수 $f'(x) = -e^{-x} - 1$ 이다.\n",
"\n",
"Newton-Raphson 식은 아래와 같다.\n",
"\n",
"$$\n",
"x_{i+1} = x_{i} - \\frac{f(x_i)}{f'(x_i)} = x_i - \\frac{e^{-x} - x}{-e^{-x} - 1}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "1b2a541b",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from matplotlib import pyplot as plt\n",
"\n",
"import numpy as np\n",
"\n",
"plt.style.use('ggplot')\n",
"plt.rcParams['figure.dpi'] = 150"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "8f6bab1a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1wAAAKPCAYAAAB5OMvxAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjErZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvzRIYmAAAAAlwSFlzAAAXEgAAFxIBZ5/SUgAAb5hJREFUeJzt3Xl81NW9//H3mUz2BMISSAJCWAyLLAoouFRRClTBavEqXUDt1fpzqVcv9XptsYpXb6WtRW0vWCt1adW2VtEqWFpRxBXZLIsiAQJhSYCEACF7JnN+f0wySUgmkJDJdzLzej4eefA957vkM/WQ8s75fs/XWGutAAAAAADtzuV0AQAAAAAQrghcAAAAABAkBC4AAAAACBICFwAAAAAECYELAAAAAIKEwAUAAAAAQULgAgAAAIAgIXABAAAAQJAQuAAAAAAgSAhcAAAAABAkBC4AAAAACBICFwAAAAAECYELAAAAAIKEwAUAAAAAQeJ2uoD2kpOTo02bNmnHjh3avn27jhw5oujoaL300kstnlddXa2///3v+uSTT5Sfny+v16vu3btr6NChmjlzprp3797knG3btmnJkiXKzs6Wx+NR3759NXXqVE2cODFInw4AAABAZxQ2gevVV1/VunXrWnXOsWPH9PDDD2vPnj1KSUnRyJEjJUkHDhzQypUrdemllzYJXGvWrNGCBQtkrdWwYcOUnJysLVu2aNGiRcrNzdUNN9zQbp8JAAAAQOcWNoErKytLmZmZGjRokAYNGqRbbrmlxeO9Xq9+8YtfaM+ePZoxY4auvfZaRUVF+fcfPHhQ8fHxjc4pKSnRokWL5PV69aMf/Ujjx4+XJB09elQPPPCAli1bprFjx2rEiBHt/wEBAAAAdDphE7iuvvrqVh3//vvva/v27Ro/fry+/e1vN9nfu3fvJn3vvfeeysrKNG7cOH/YkqSUlBTNmjVLjz32mJYuXUrgAgAAACApghfNWLFihSTpyiuvPOVz1q9fL0maMGFCk31jxoxRdHS0Nm/erKqqqvYpEgAAAECnFjYzXK1RXl6unTt3Kj4+XoMHD1Z2drbWrVunkpIS9ezZU+PGjVO/fv2anLdnzx5J0sCBA5vsc7vd6tevn3bu3Km8vDxlZmaedp0HDhyQtfa0r4PT17NnT0lSYWGhw5Ug1DA20BLGBwJhbCAQxkboMcYoLS2tzedHZODat2+frLVKS0vTc889p3/84x+N9v/lL3/RlVdeqVmzZvn7ysrKVFpaKknNrlxY179z504VFha2S+Cy1hK4Qgz/PRAIYwMtYXwgEMYGAmFshI+IDFx1wWnPnj3atWuXrrzySk2dOlVxcXFau3atnn/+eb355pvq1auXpkyZIkmqqKjwnx8bG9vsdev6Gx7bkjlz5jTpi4mJ0fz58yXV/4YDznO7fX9VUlNTHa4EoYaxgZYwPhAIYwOBMDbCT0Q+w+X1eiVJNTU1uvDCCzV79mz16tVLXbp00aRJk/S9731PkvT66687WSYAAACATi4iZ7ji4uL825deemmT/Zdeeqmee+45HT58WAcOHFBaWlqjcyorK5WQkNDkvMrKyibXb8mCBQta3F9YWMh0coio+y1TQUGBw5Ug1DA20BLGBwJhbCAQxkboMcYoPT29zedH5AxXr169/NvNTdfGxsaqS5cuknwvR5akhIQEf8gqKipq9rp1/dwKCAAAAECK0MDVs2dPJScnS/K9zPhEXq/X/5xXw9mq/v37S5JycnKanOPxeLRnzx5FR0crIyMjGGUDAAAA6GQiMnBJ0tixYyVJX3zxRZN92dnZ8ng8iomJUZ8+ffz9Y8aMkSStXr26yTkbNmxQdXW1RowYoZiYmCBVDQAAAKAzidjA9c1vflMul0tvvvmmdu3a5e8/duyYnnvuOUm+Z7nqVoqRpEmTJik+Pl7r1q3TZ5991uicF198UZI0ffr0DvoEAAAAkavu9Tl88dXar45mrBPfNQg2bNig1157zd/evn27jDEaPHiwv++aa67xz1JJ0t///nc999xzio6OVlZWlmJjY7Vt2zaVlpZqwIABmjdvnuLj4xt9n9WrV+vxxx+XJA0fPlzJycnavHmzSktLdfnll+v73/9+u32m/Px8RwYFmuIBVgTC2EBLGB8IhLHRNh6PR+Xl5aqqqgrbfyO5XL75kLpVtdH+jDGKiYlRfHx8o8mVlo4/nUUzwmaVwuLiYm3fvr1Rn7W2UV9xcXGj/ZdffrkyMjL01ltvaceOHaqurlbv3r01bdo0XXnllc2+b2vChAl66KGHtGTJEm3fvl0ej0d9+vTR1KlTm13xEAAAAKfP4/Ho2LFjiouLU0pKij+YhJu6AODxeByuJHx5vV5VVFTo2LFj6tq16ymFrtMRNjNc4YgZrtDBbyIRCGMDLWF8IBDGRusdP35cLpdLiYmJTpcSVASujlNaWiqv1+tfTC8QloUHAABA2Kuqqjrld50CpyI2NlZVVVVB/z4ELgAAAIS0usUOwvU2QjgjKiqqQxbSYNQCAAAAQJAQuAAAAAAgSAhcAAAAABAkBC4AAAAAAX3wwQfq06eP1q5d6+/75z//qUGDBumdd95xsLLOgcAFAAAAIKCLL75YEydO1KOPPipJ+vTTT3X77bfrl7/8pSZPnuxwdaEvbF58DAAAAKBlR44c0dGjR1s8Jjk5WT179mzUd//992vKlCn69a9/rUWLFun+++/XjBkzglhp+CBwAQAAABHi2Wef1YIFC1o85tprr9UTTzzRqG/YsGG66qqr9POf/1z33HOPbrzxxuAVGWa4pRAAAADo5L766ivdc889mjBhggYMGKDRo0frrrvuUn5+fqPj/vM//1O5ubktfjUXyL744gu99957io6OVo8ePTrqYzVxxx13qE+fPnryySeb7FuzZo0GDhyoUaNGaffu3R1fXAAELjRhDx+S9y+/V81vHlbN4w84XQ4AAABa8Pzzz2vq1Kn6y1/+ol69emny5MlKSkrSq6++qm9+85sqKiryH+tyueR2u1v8OvEF07t27dL3vvc9XX/99frhD3+oBQsWqLS0tKM/piTpnnvukdvt1u9+9zsdP37c379jxw59//vfl8vl0h/+8AdlZmY6Ul9zCFxoqrpKdsXfpE1rpa2bZD3VTlcEAACAZvz973/X/fffr379+mn58uV688039bvf/U4ffPCBvvOd7ygvL09PPfVUm6+fn5+v73znO5o6daruu+8+3XrrraqpqTmta56OAQMGaObMmTp69KgWL14sSSooKNDs2bN1/Phx/fa3v9XZZ5/tSG2BELjQVM/ekqkdGtYrFR50th4AAAA0UVpaqvvuu09JSUl66aWXdNZZZ/n3RUVFae7cuZKkDz/8sE3XLyoq0ne/+12NGjXKv0JhUlKS/uM//kNPP/20Dh505t+Id999t+Li4vTMM8/owIEDuuGGG7Rnzx7Nnz9fX//61x2pqSUsmoEmjDta6tlLKjjg6ziYL6X1dbYoAACAZlhrpXJnbm87LfGJMsac1iVefPFFFRYW6o477lC/fv2a7O/WrZuSkpJOuiphIN27d9fKlSub9P/gBz/QD37wg5Oe/4Mf/EDbtm1r1fd88skndc4557R4TEZGhmbPnq1nnnlGkydPVlFRkebMmaPvfve7rfpeHYXAheb1zvAHLntwv4zOdbggAACAZpSXyntXaP5DuyWuJ1+WEpJO6xr/+Mc/JEkLFy7UwoULAx7Xv3//0/o+bbVv3z7t3LmzVeeUl5ef0nG33HKLFi9erKKiIl133XX60Y9+1JYSOwSBC80yvTJktcHXOJTnbDEAAABo4ssvv1RcXJyuvPLKFo9reKthR/r73/8elOtaa/XQQw/5Zjclud2hHWlCuzo4p3eGf9Meym/hQAAAAHS0qqoqHT9+XP369WvyzqxwN2/ePC1dulSTJ0/W559/rldeeUW33XabBg4c6HRpzSJwoVm+Ga5aB/c7WQoAAEBg8Ym+2/M6m/jE0zo9JiZGcXFxysvLU1lZmRISEtqpsND29NNPa/HixTrnnHP01FNP6Q9/+IP+53/+R4899pgWLVrkdHnNYpVCNK/BDJeKCmWrKp2rBQAAIABjjExCUuf7Os0FMyTp4osvlsfj0X333dfk2SdrrT788EOtXbv2tL9PqHjzzTf18MMPKzMzUy+88ILi4+N1/fXXq1evXnrzzTf15ZdfOl1iswhcaF6PVCmqwQRo3YqFAAAACAlz585VSkqKXnvtNY0fP16zZ8/WrbfeqquvvlrnnHOOvv3tbzv2guL29umnn+ruu+9W9+7d9dJLL6lHjx6SpPj4eN1xxx2y1urnP/+5w1U2j8CFZhlXlJSaVt/BbYUAAAAhZfDgwfrnP/+pWbNmKTExUR999JHee+89HTp0SGPHjtUvf/lLjR8/3ukyT1t2drZuuukmuVwuvfDCC8rMzGy0f9asWUpLS9OKFSu0bt06Z4psAc9wIbDeGdKBfZIkezBfpz/xDQAAgPbUp0+fkJ3ZaS9ZWVkt3i4YFxen9evXd2BFrcMMFwIyvdLrGywNDwAAALQagQuB9e7j37TcUggAAAC0GoELATWe4eJdXAAAAEBrEbgQWMOl4Y8dka0oc64WAAAAoBMicCGwlB5STEx9+yCzXAAAAEBrELgQkHG5pNT62wotC2cAAAAArULgQssa3lZ4kMAFAAAAtAaBCy0yBC4AAACgzQhcaFmv+sDFLYUAAABA6xC40CLTIHDx8mMAAOAEY4wkyVrrcCUIJ3XjqW58BQuBCy1LaxC4So7Llh53rhYAABCx3G63qqurnS4DYaS6ulputzvo34fAhZYlp0hx8fVtnuMCAAAOiImJUUVFBbNcaBfWWlVUVCim4SuQgiT4kQ6dmjHG9xzXnp2SJHswT2bgEIerAgAAkSY+Pl6VlZUqKSlRfHy8oqKinC4pKOoCJcEyeGpqalReXi6v16v4+PiTn3CaCFw4KdM7Q7Y2cOngfmeLAQAAEckYo65du6q8vFzHjh0L20DicvluQPN6vQ5XEr6MMYqLi1PXrl2D/vyWRODCqUjrU799gMAFAACc4XK5lJiYqMTExLANXKmpqZKkgoIChysJXx0RshoicOHketcHLssMFwAACAEd/Y/mjlL3ucL180UiFs3ASZm0vvWNg3myTHEDAAAAp4TAhZPr3WBp+OoqqYgpbgAAAOBUELhwUiYuXurWs76D57gAAACAU0LgwqlJ4zkuAAAAoLUIXDglpnfDlQr3OVcIAAAA0IkQuHBqGiycYbmlEAAAADglBC6cEsO7uAAAAIBWI3Dh1DQMXEcPy1aUOVcLAAAA0EkQuHBquvWUYmLq2wfznKsFAAAA6CTcThfQXnJycrRp0ybt2LFD27dv15EjRxQdHa2XXnrplK/x8MMPa/PmzZKk3/3ud0pJSWn2uG3btmnJkiXKzs6Wx+NR3759NXXqVE2cOLEdPkloMi6X1KuPtG+XJN9zXKb/YIerAgAAAEJb2ASuV199VevWrWvz+e+//742b94sY4ystQGPW7NmjRYsWCBrrYYNG6bk5GRt2bJFixYtUm5urm644YY21xDqTFof2drAxXNcAAAAwMmFTeDKyspSZmamBg0apEGDBumWW2455XOLi4v1xz/+UaNHj1ZeXp4KCgqaPa6kpESLFi2S1+vVj370I40fP16SdPToUT3wwANatmyZxo4dqxEjRrTLZwo5aSwNDwAAALRG2DzDdfXVV+u6667T2LFjA94KGMjzzz+viooK3XTTTS0e995776msrEzjxo3zhy1JSklJ0axZsyRJS5cubXXtnUaDd3GxNDwAAABwcmETuNrqX//6lz766CPNmDFDaWlpLR67fv16SdKECROa7BszZoyio6O1efNmVVVVBaVWpzVaGv7Qflmv17liAAAAgE4gogNXZWWlnnnmGfXp00dXXXXVSY/fs2ePJGngwIFN9rndbvXr10/V1dXKywvTFfwaBq6qKunIYedqAQAAADqBiA5cf/nLX1RQUKCbb75ZbnfLj7OVlZWptLRUktS9e/dmj6nrLywsbN9CQ4SJS5BSGnx2nuMCAAAAWhQ2i2a0Vk5Ojt5++21dcsklOuuss056fEVFhX87Nja22WPq+hse25I5c+Y06YuJidH8+fMlST179jyl63SkojMGqOpokSQpseSYElNTHa6oY9QF8tQI+bw4dYwNtITxgUAYGwiEsRF+InKGy+v16umnn1ZiYqJmz57tdDmdSlSffv7tmrxcBysBAAAAQl9EznAtW7ZMu3bt0q233qouXbqc0jlxcXH+7crKSiUkJDQ5prKyssmxLVmwYEGL+wsLC1t8J5gTvF17+LfLdu1QZYAl9MNN3W+ZAr0yAJGLsYGWMD4QCGMDgTA2Qo8xRunp6W0+PyID1/r162WM0apVq/TBBx802nf06FFJ0mOPPSa3261vf/vbGjp0qBISEpSQkKCysjIVFRU1G7iKiny32oXirYDtxaT1kT8C5vMMFwAAANCSiAxckmSt1datWwPuz87OluR7KXKd/v37a+vWrcrJyVHfvn0bHe/xeLRnzx5FR0crIyMjOEWHgvQz6rePHpYtL5OJbxo+AQAAAERo4Jo3b17AfXfccYcKCgr0u9/9rskLlMeMGaOtW7dq9erVuvjiixvt27Bhg6qrq3XOOecoJiYmCFWHiG49pdg4qbJ2YZAD+6QBWc7WBAAAAISoiFw0o60mTZqk+Ph4rVu3Tp999pm//9ixY3rxxRclSdOnT3eqvA5hXC4prX52z+btdbAaAAAAILSFzQzXhg0b9NprrzXq83g8mjt3rr99zTXXaMyYMW3+HklJSbrtttv0+OOPa8GCBRo+fLiSk5O1efNmlZaW6vLLL9fIkSPbfP3OwqT3lc3d4WvkE7gAAACAQMImcBUXF2v79u2N+qy1jfoaPo/VVhMmTNBDDz2kJUuWaPv27fJ4POrTp4+mTp2qSy+99LSv3yk0nOHi5ccAAABAQGETuCZOnKiJEyee9nUWLlx40mOGDh2qn/zkJ6f9vTork9GvwUqFzHABAAAAgfAMF1ovvcEKjQUHZaurnKsFAAAACGEELrRearoUVTs5ar3Swf3O1gMAAACEKAIXWs1ERUm96t+2bXkBMgAAANAsAhfaJqPBC5B5jgsAAABoFoELbWLSGwYuZrgAAACA5hC40DYNl4ZnhgsAAABoFoELbWIy+tU3Du6XralxrhgAAAAgRBG40Da9MyRjfNsej1R40Nl6AAAAgBBE4EKbmJhYqWfv+g5uKwQAAACaIHCh7Ro9x8XCGQAAAMCJCFxoM9Noafg9zhUCAAAAhCgCF9quwdLwzHABAAAATRG40GamwS2FOrBP1lrnigEAAABCEIELbdfw5ccV5dKRQudqAQAAAEIQgQttZhISpZTu9R15rFQIAAAANETgwulp8AJkm8fCGQAAAEBDBC6cFpPRv75B4AIAAAAaIXDh9DRYGp4ZLgAAAKAxAhdOi2lwS6Hy9rJSIQAAANAAgQunp2HgqiyXigqcqwUAAAAIMQQunBYTnyB1T63v2J/rXDEAAABAiCFw4fSxUiEAAADQLAIXTlvj57gIXAAAAEAdAhdOX58GM1z7CVwAAABAHQIXTlujGa4De2W9XueKAQAAAEIIgQunL73+XVyqqpIKDzpXCwAAABBCCFw4bSY2TkpNq+/IY6VCAAAAQCJwob1k8BwXAAAAcCICF9qFyWhwW2HeXucKAQAAAEIIgQvtI6O/f9NySyEAAAAgicCFdtJ4pcJ9sjU1zhUDAAAAhAgCF9pHel/J1A4nj0cqyHe2HgAAACAEELjQLkx0jNQrvb4jj4UzAAAAAAIX2k+DhTNYqRAAAAAgcKEdmT71C2doPwtnAAAAAAQutJuGgcvu3+1cIQAAAECIIHCh/fTNrN8+mC9bVelYKQAAAEAoIHCh/fRKl6JjfNvWK+XzAmQAAABENgIX2o1xRUkN3sdl9+12rhgAAAAgBBC40K4aLZyxj4UzAAAAENkIXGhfDZ7jYuEMAAAARDoCF9pV4xmu3Y7VAQAAAIQCAhfaV8OVCo8fky0+4lgpAAAAgNMIXGhXpkuKlNy1voPnuAAAABDBCFxofw2f4+K2QgAAAEQwAhfanemTWd/YzwwXAAAAIpfb6QLaS05OjjZt2qQdO3Zo+/btOnLkiKKjo/XSSy81Odbr9Wrbtm1av369vvzySx06dEhlZWXq0aOHRo4cqauvvlq9evUK+L22bdumJUuWKDs7Wx6PR3379tXUqVM1ceLEIH7CToQZLgAAAEBSGAWuV199VevWrTulYw8dOqQHH3xQktS9e3dlZWXJ5XJpx44dWrFihT7++GP9+Mc/1tChQ5ucu2bNGi1YsEDWWg0bNkzJycnasmWLFi1apNzcXN1www3t+rk6I9O3v2xdI3+vbE2NTFSUkyUBAAAAjgibwJWVlaXMzEwNGjRIgwYN0i233NLi8aNHj9a3vvUtDR8+3N9XXV2tZ555Ru+//75+/etf69e//rXc7vr/iUpKSrRo0SJ5vV796Ec/0vjx4yVJR48e1QMPPKBly5Zp7NixGjFiRHA+ZGeRfoZkXJL1StVV0qF8Kb2v01UBAAAAHS5snuG6+uqrdd1112ns2LFKSUlp8di0tDTNnTu3UdiSpOjoaN18881KSEhQYWGhsrOzG+1/7733VFZWpnHjxvnDliSlpKRo1qxZkqSlS5e2zwfqxExMrNQ7vb6DFyADAAAgQoVN4GovMTExSk/3hYWioqJG+9avXy9JmjBhQpPzxowZo+joaG3evFlVVVXBLzTENVw4g+e4AAAAEKkIXCfwer0qLCyUpCYzZXv27JEkDRw4sMl5brdb/fr1U3V1tfLy8oJeZ8jr29+/aVmpEAAAABGKwHWCjz/+WMeOHVOXLl00ZMgQf39ZWZlKS0sl+RbaaE5df11gi2SmwUqFYoYLAAAAESpsFs1oD4WFhXr++eclSdddd52io6P9+yoqKvzbsbGxzZ5f19/w2JbMmTOnSV9MTIzmz58vSerZs+cpXScUeUaNlT92Fh5Uj8QEuRISnSzptNQtnpKamupwJQg1jA20hPGBQBgbCISxEX6Y4apVUVGhxx57TMePH9e5556rKVOmOF1SpxbVK12mQcDy7N7hYDUAAACAM5jhkuTxePSrX/1KOTk5Gjp0qO66664mx8TFxfm3KysrlZCQ0OSYysrKJse2ZMGCBS3uLywslLW2xWNCme3TX9r+pSTpyOYNcqVmOFxR29X9lqmgoMDhShBqGBtoCeMDgTA2EAhjI/QYY/yL6rVFxM9web1e/eY3v9HGjRvVv39//fd//7diYmKaHJeQkOAPWSeuXlinrr8z3wrYnswZDRYX2bvLuUIAAAAAh0R84Fq8eLE+/fRTpaen6/7771diYuDnjPr39628l5OT02Sfx+PRnj17FB0drYyMzjuT067OGODftAQuAAAARKCIDlwvv/yyVqxYoZ49e+qnP/2punbt2uLxY8aMkSStXr26yb4NGzaourpaI0aMaHaGLBI1muHanyvr8ThXDAAAAOCAiA1cS5cu1RtvvKGUlBT99Kc/PaXbACdNmqT4+HitW7dOn332mb//2LFjevHFFyVJ06dPD1rNnU7GGVJUlG/bUy0d3O9sPQAAAEAHC5tFMzZs2KDXXnutUZ/H49HcuXP97WuuuUZjxozR7t279cc//lGS1KtXLy1ZsqTZa06aNElDhw71t5OSknTbbbfp8ccf14IFCzR8+HAlJydr8+bNKi0t1eWXX66RI0cG4dN1TiY6RkrrK9W++NjuzZHp0/8kZwEAAADhI2wCV3FxsbZv396oz1rbqK+4uFiSVFpa6l/9Lzs7W9nZ2c1e86yzzmoUuCRpwoQJeuihh7RkyRJt375dHo9Hffr00dSpU3XppZe250cKC+aMgbK1gUt7d0sTHC0HAAAA6FDGduZ1x8Ncfn5+p14WXpK8/3xD9q/P+hrDRitqzsPOFtRGLNGKQBgbaAnjA4EwNhAIYyP0sCw8QpppsFKh9u7q9AESAAAAaA0CF4Krb4PAVVIsHW3+HWYAAABAOCJwIahMchepW4MVIPc2fYcZAAAAEK4IXAg+XoAMAACACEXgQtCd+BwXAAAAECkIXAg6c8ZA/zYzXAAAAIgkBC4EX8MZroJ82Yoy52oBAAAAOhCBC8HXs7cUF+/btlbat9vRcgAAAICOQuBC0BmXq9Hy8NxWCAAAgEhB4EKHMP0H1TdydzhXCAAAANCBCFzoGP0aLJyRy7u4AAAAEBkIXOgQpl+DGa78PbLVVc4VAwAAAHQQAhc6RvoZUnSMb7umRtqf62w9AAAAQAcgcKFDmKgoqW+mv21zdzpXDAAAANBBCFzoMI0WzthD4AIAAED4I3Ch4zR4josZLgAAAEQCAhc6TKOFM/bvlvV4nCsGAAAA6AAELnScPv2kKLdv2+OR8vY4Ww8AAAAQZAQudBjjjpb69Pe3Lc9xAQAAIMwRuNChWDgDAAAAkYTAhY7Vb6B/k4UzAAAAEO4IXOhQpv/g+sa+XbI1Nc4VAwAAAAQZgQsdq09/yVU77KqqpAP7na0HAAAACCICFzqUiYmV0s/wt1k4AwAAAOGMwIUO1+h9XLk7nCsEAAAACDICFzpeg+e4WDgDAAAA4YzAhQ5nMhssnLFnJwtnAAAAIGwRuNDx+g5osHBGpZS/19l6AAAAgCAhcKHDmdhYKaO/v213b3ewGgAAACB4CFxwhBlwZn2DhTMAAAAQpghccEaD57jsLma4AAAAEJ4IXHCEyWwww7Vvt2x1tXPFAAAAAEFC4IIzMvpL0TG+7RqPtG+3o+UAAAAAwUDggiOM2y2dMcDfZuEMAAAAhCMCFxxjGrwAWQQuAAAAhCECF5zT4Dkuy0qFAAAACEMELjim0dLweXtlKyucKwYAAAAIAgIXnNO7jxQX79u2Xil3p7P1AAAAAO2MwAXHGJdLavAcFwtnAAAAINwQuOAok8nCGQAAAAhfBC44quELkJnhAgAAQLghcMFZDZeGLzggW3rcuVoAAACAdkbggrN69paSu9a3d2U7VwsAAADQzghccJQxRhqQ5W/bHAIXAAAAwgeBC44zDQPXrm0OVgIAAAC0LwIXHGcGDqlv5GTLWutcMQAAAEA7InDBeZlnSsb4tstKpIN5ztYDAAAAtBO30wW0l5ycHG3atEk7duzQ9u3bdeTIEUVHR+ull15q8bxVq1Zp+fLl2rdvn9xut7KysjRjxgwNGTIk4Dnbtm3TkiVLlJ2dLY/Ho759+2rq1KmaOHFiO3+qyGASEqW0vlL+XkmS3ZUtk9bH4aoAAACA0xc2gevVV1/VunXrWnXOCy+8oGXLlikmJkajRo1SdXW1Nm3apI0bN2rOnDk677zzmpyzZs0aLViwQNZaDRs2TMnJydqyZYsWLVqk3Nxc3XDDDe31kSKKGZglWxu4lLNNOv9SZwsCAAAA2kHYBK6srCxlZmZq0KBBGjRokG655ZYWj9+yZYuWLVum5ORkPfLII0pPT5ckZWdna968eVq0aJGGDx+upKQk/zklJSVatGiRvF6vfvSjH2n8+PGSpKNHj+qBBx7QsmXLNHbsWI0YMSJ4HzRcDRwiffyuJN8MFwAAABAOwuYZrquvvlrXXXedxo4dq5SUlJMe/9Zbb0mSZsyY4Q9bki+4TZ48WWVlZVq5cmWjc9577z2VlZVp3Lhx/rAlSSkpKZo1a5YkaenSpe3waSKPGdDgFs59u2SrKp0rBgAAAGgnYRO4WqOqqkpbtmyRJE2YMKHJ/rq+9evXN+qvazd3zpgxYxQdHa3NmzerqqqqvUsOfxn9pJhY33ZNjbRnp7P1AAAAAO0gIgNXXl6eqqur1aVLF/Xo0aPJ/gEDBkiScnNzG/Xv2bNHkjRw4MAm57jdbvXr10/V1dXKy2OVvdYyUVG+1Qpr2RzexwUAAIDOLyIDV2FhoSQ1G7YkKS4uTomJiSotLVV5ebkkqaysTKWlpZKk7t27N3teXX/d9dE6DV+ArBye4wIAAEDnFzaLZrRGRUWFJCkmJibgMbGxsSotLVVFRYXi4+P959TtC3ROw+ufzJw5c5r0xcTEaP78+ZKknj17ntJ1wkXFOefq6D+WSJJM7g6lpqY6XFE9t9v3VyWUakJoYGygJYwPBMLYQCCMjfATkTNc1lpJkql72W4Lx6DjRGed5d/2Fh5UTVGBg9UAAAAApy8iZ7ji4+MlSZWVgVfCq1v4Ii4urtGfdeclJCQ0Oafueg2PbcmCBQta3F9YWBh5wa97T6nId0vm4bWfyow53+GCfOp+y1RQQAhEY4wNtITxgUAYGwiEsRF6jDGNVjVvrYic4aq7Ve/w4cPN7q+oqFBpaakSExP94SwhIcEfsoqKipo9r64/0m4FbFcNnuOyOV85WAgAAABw+iIycGVkZCg6OlrFxcXNhq5du3ZJkvr169eov3///pKknJycJud4PB7t2bNH0dHRysjICELVkcEMGubftjsJXAAAAOjcIjJwxcTEaMSIEZKk1atXN9lf1zd27NhG/WPGjAl4zoYNG1RdXa0RI0a0uBgHWmYGDa1v7N4hW13tXDEAAADAaYrIwCVJ06ZNkyQtWbJE+fn5/v7s7GytWLFC8fHxuuyyyxqdM2nSJMXHx2vdunX67LPP/P3Hjh3Tiy++KEmaPn16B1QfxvoNlKJrA6unmhcgAwAAoFMLm0UzNmzYoNdee61Rn8fj0dy5c/3ta665xj9LNWrUKF1xxRV6++23de+992rkyJGqqanRpk2b5PV6deeddyopKanR9ZKSknTbbbfp8ccf14IFCzR8+HAlJydr8+bNKi0t1eWXX66RI0cG/8OGMeOOljIHS9u/lCTZnVsbz3oBAAAAnUjYBK7i4mJt3769UZ+1tlFfcXFxo/033nijMjMztXz5cm3evFlRUVEaMWKErrnmGg0d2vw/8idMmKCHHnpIS5Ys0fbt2+XxeNSnTx9NnTpVl156aft/sAhkBg2T9QcunuMCAABA52VsxK073nnk5+dH3rLwkuzGNfL+3yO+RpcUuR57ocV3pnUElmhFIIwNtITxgUAYGwiEsRF6WBYe4Wdgg9nF4qNS4UHHSgEAAABOB4ELIcckd5F69/G37c6tDlYDAAAAtB2BCyHJDG4wy7WDwAUAAIDOicCF0MQLkAEAABAGCFwISY2Wgt+fK1te5lwxAAAAQBsRuBCa0vpKCbXvQbNWytnmbD0AAABAGxC4EJKMyyU1mOVi4QwAAAB0RgQuhKyGtxVaFs4AAABAJ0TgQsgyg+sXzlBOtmxNjXPFAAAAAG1A4ELoysySoty+7cpyaW+Os/UAAAAArUTgQsgysbFS5mB/227/0sFqAAAAgNYjcCGkmcHD/dt2+xcOVgIAAAC0HoELIc2ceVZ9Y/uXstY6VwwAAADQSgQuhLbBwyRjfNslxdKBfc7WAwAAALQCgQshzSQmSX36+9vcVggAAIDOhMCFkGfOrH+OS9kELgAAAHQeBC6EvjNH+DdZqRAAAACdCYELIa/RDFdRgezhQ84VAwAAALQCgQshz6R0l1LT/G2e4wIAAEBnQeBCp2CyGi8PDwAAAHQGBC50Dg3ex8VzXAAAAOgsCFzoFBo9x5W/V/b4MeeKAQAAAE4RgQudQ2q61LV7fZvl4QEAANAJELjQKRhjGj3HZbO3OFgNAAAAcGoIXOg8ho70b9ptmx0sBAAAADg1BC50GiarPnBpfy7PcQEAACDkEbjQefTO4DkuAAAAdCoELnQaxhiZIQ1vK9zkYDUAAADAyRG40LkMGeHftNtYOAMAAAChjcCFTsU0WDhDeXtki486VgsAAABwMgQudC6p6VJKj/o2y8MDAAAghBG40Kn4nuPitkIAAAB0DgQudD5DeB8XAAAAOgcCFzqdhisVKn+vbPER54oBAAAAWkDgQueTmiZ17+lv2m28jwsAAAChicCFTscYI5PVYJaL93EBAAAgRBG40Dk1WB7ebiVwAQAAIDQRuNApmaGj6xuH8mQPH3KuGAAAACAAAhc6JdMjVerdx9+2Wzc6WA0AAADQPAIXOi0zbFR9g9sKAQAAEIIIXOi0zLD62wrt1n/JWutgNQAAAEBTBC50XkNGScb4to8fk/bnOlsPAAAAcAICFzotk5gk9Rvkb9uveI4LAAAAoYXAhU7NDG9wW+GXBC4AAACEFgIXOrVGy8Nnb5H1eJwrBgAAADgBgQud2+Bhkjvat11ZIe3KdrYeAAAAoAECFzo1ExMrnTnc3+Z9XAAAAAglbqcLcFp2drbefPNNbdu2TSUlJYqLi9OAAQM0ZcoUTZgwodlzVq1apeXLl2vfvn1yu93KysrSjBkzNGTIkA6uHpJkho7yBy27daP0ze84XBEAAADgE9GB69NPP9UTTzwha60GDRqks846S0eOHNEXX3yhLVu26KqrrtL3vve9Rue88MILWrZsmWJiYjRq1ChVV1dr06ZN2rhxo+bMmaPzzjvPoU8Tucyws2Vf/6OvsWubbHmZTHyCs0UBAAAAiuDAVVNTo9///vey1uruu+/WBRdc4N+XnZ2thx56SG+++aYmTZqktLQ0SdKWLVu0bNkyJScn65FHHlF6err/+Hnz5mnRokUaPny4kpKSHPlMEav/QCkxWSo9LtXUSNs2SWc3PzsJAAAAdKSIfYZr//79Ki4uVp8+fRqFLUnKysrS6NGjZa1VTk6Ov/+tt96SJM2YMcMftuqOnzx5ssrKyrRy5cqO+QDwM64omWENlof/4nMHqwEAAADqRWzgio6OPqXj6marqqqqtGXLFklq9tmuur7169e3U4VolbPO8W8SuAAAABAqIjZw9e7dW71799b+/fv1ySefNNqXnZ2tjRs3qlevXho+3LcCXl5enqqrq9WlSxf16NGjyfUGDBggScrNzQ1+8WjCDK8PXCo4IHsoz7liAAAAgFoR+wyXy+XS7bffrp///Od64okn9NZbb6l37946cuSIvvrqKw0ePFh33nmn3G7f/0SFhYWS1GzYkqS4uDglJiaqtLRU5eXlio+PP2kNc+bMadIXExOj+fPnS5J69uzZ1o8XeVJTVdhvoDx7fLeAJuZuV+JZo09y0qmrGwepqantdk2EB8YGWsL4QCCMDQTC2Ag/ETvDJUnDhg3TvHnz1KtXL+3cuVOffPKJtm7dqri4OI0cOVLdunXzH1tRUSHJF4gCiY2NbXQsOlbM2eP921Wff+ZgJQAAAIBPxM5wSdJHH32kp556Smeeeabuvvtu9e3bV0eOHNFbb72lJUuWaMuWLZo3b57cbrestZIkY0zA69Udc6oWLFjQ4v7CwsJWXzOS2YFD/duVm9brUH6ejPvUntU7mbrfMhUUFLTL9RA+GBtoCeMDgTA2EAhjI/QYYxotmNdaETvDlZ+fr4ULF6pLly667777NHjwYMXFxSk9PV233HKLxo4dq+zsbL3//vuS5L9FsLKyMuA1q6qqJPluL4QDzhwuRdfOQFaWSzu/crYeAAAARLyIDVwff/yxampqNHr06GYD0vnnny9J+uKLLyTVP091+PDhZq9XUVGh0tJSJSYmntLzW2h/JiZWyjrL37ZfbHCwGgAAACCCA1dRUZEkKSEhodn9df0lJSWSpIyMDEVHR6u4uLjZ0LVr1y5JUr9+/YJRLk6ROWuMf5vl4QEAAOC0iA1cKSkpkqSdO3c2u3/Hjh2S6u+jjYmJ0YgRIyRJq1evbnJ8Xd/YsWPbu1S0gmnwPi7tyZEtPuJcMQAAAIh4ERu4xo0bJ0naunWr/vnPfzbal52drWXLlklq/JLjadOmSZKWLFmi/Pz8RsevWLFC8fHxuuyyy4JdOlqSfobUrX45ffvFv5yrBQAAABEvYlcpHDhwoK688kq99dZbWrx4sf7xj3+oT58+OnLkiLKzs2Wt1de//nWNGjXKf86oUaN0xRVX6O2339a9996rkSNHqqamRps2bZLX69Wdd96ppKQkBz8VjDEyI8bIflgbojevk86/1NmiAAAAELEiNnBJ0uzZszVkyBC98847ysnJUV5enuLi4jR8+HBNmjRJF110UZNzbrzxRmVmZmr58uXavHmzoqKiNGLECF1zzTUaOnRoM98FHc2MHOcPXPaLDbI1NTJRUQ5XBQAAgEgU0YFLks477zydd955rTpn4sSJmjhxYnAKwukbNlpyuyWPRyor9S0P32D1QgAAAKCjROwzXAhfJi5eyhrhb9st6xysBgAAAJGMwIWwZEbWrxZpNxG4AAAA4AwCF8KSGXlufWN/ruzhAueKAQAAQMRq98B155136o033tCxY8fa+9LAKTO9M6Re6f623cwsFwAAADpeuweuQ4cO6U9/+pNuu+02LViwQJs2bWrvbwGcEjNynH/bblnvYCUAAACIVO0euL71rW+pe/fuqqmp0Weffab//d//1Z133qm//e1vzHqhQ5lR9YFLWzfKVlc5VwwAAAAikrHW2va+qNfr1eeff6533nlHGzdulNfrlSRFRUXp3HPP1aRJkxq9UBjNy8/PVxD+80QMW10t739+T6qskCS57npQZsTYk5zVvNTUVElSQQHPgqExxgZawvhAIIwNBMLYCD3GGKWnp5/8wACC8h4ul8ulsWPHauzYsSoqKtJ7772nlStXqrCwUKtXr9bq1avVq1cvff3rX9fEiRPVtWvXYJSBCGeio33v5PrXZ5J8qxW2NXABAAAAbRGUGa7mWGu1ceNGrVixQhs2bFBNTY0kZr1awgzX6fN+8A/ZPy70NXr0kuvRZ2SMafV1+G0TAmFsoCWMDwTC2EAgjI3QE5IzXM0xxujss8/W2WefraKiIv3617/W1q1bVVNT45/1SktL05VXXqnLLrtMLhcr1uP0mVHj5I+shw9J+3dLfQc4WBEAAAAiSYemmsLCQr3yyiuaO3eutm7d6u/PzMyUy+XSgQMH9Mwzz2ju3LkqLi7uyNIQpkxKDynzTH/b/muNg9UAAAAg0gR9hsvr9WrdunV69913tWnTJv8CGklJSZo4caImT56stLQ0HT16VP/85z+1bNky5eTk6OWXX9att94a7PIQAczZ42V3b5ck2Y1rpOkzHa4IAAAAkSJogevQoUN699139f777+vo0aP+/qysLE2ZMkXnn3++3O76b5+SkqLrrrtOY8eO1U9+8hN9/vnnwSoNEcaMPk/2jRd9jd3bZY8e9s18AQAAAEHW7oFr9erVWrFihbZs2eJf8CE+Pl4XXXSRpkyZon79+rV4/qBBg5SSktIopAGnpU9/qWdvqfCgJMluXCtzyTccLgoAAACRoN0D1+OPP+7fzszM1OTJk3XRRRcpLi7u1Ityd9haHogAxhjfLNe7b0mS7L8+kwhcAAAA6ADtnmyio6N1wQUXaMqUKRo8eHCbrrFw4cJ2rgqRrmHg0lcbZSvKZeLinS0KAAAAYa/dA9fTTz+txMTE9r4scHrOPEtKSJTKSiWPR/ric2nsBU5XBQAAgDDX7svCE7YQiozbLTNinL9tN37mYDUAAACIFLxdGJHj7PH+TbtpnWxNjYPFAAAAIBIQuBAxzIgxUlTtXbSlx6UdXzpbEAAAAMIegQsRw8QnSMNG+dt2w6cOVgMAAIBIQOBCRDHnnO/fths+lfV6HawGAAAA4Y7AhYhizh4vmdphf/SwtHu7swUBAAAgrBG4EFFMlxTpzOH+NrcVAgAAIJgIXIg4ZkzD2wo/kbXWwWoAAAAQzghciDjmnAn1jYID0v7djtUCAACA8EbgQsQx3VOlAVn+NrcVAgAAIFgIXIhIjW8rJHABAAAgOAhciEgNA5f258oe2O9cMQAAAAhbBC5EJNMrQ+qb6W/bDZ84VwwAAADCFoELEcuMucC/bdd/7GAlAAAACFcELkQsM+7C+saeHNmDec4VAwAAgLBE4ELEMulnSH36+9t23UcOVgMAAIBwROBCRDPnfs2/bdd+6GAlAAAACEcELkQ0M+6i+sb+XNn8vc4VAwAAgLBD4EJEM70zpH4D/W27ltsKAQAA0H4IXIh4ZlyD2wp5jgsAAADtiMCFiNdotcL8vbL7c50rBgAAAGGFwIWIZ1LTpMwz/W1muQAAANBeCFyAGi+eYdd8KGutg9UAAAAgXBC4AJ2wWuGhPCl3h3PFAAAAIGwQuABJpkeqdOZwf9t+tsrBagAAABAuCFxALTN+on/brv1Q1lvjXDEAAAAICwQuoJYZe4EUFeVrHDsifbXJ2YIAAADQ6RG4gFomqYs0Yqy/bVdzWyEAAABOD4ELaMCMv8S/bT//VLaq0sFqAAAA0Nm5nS4gFBw9elRvvPGGPv/8cxUWFiomJka9evXSyJEjNWvWrCbHr1q1SsuXL9e+ffvkdruVlZWlGTNmaMiQIQ5Uj/ZkRp0nGxsvVZZLFeWyG9fKnHvRyU8EAAAAmhHxM1zZ2dn6z//8T7399tuKiorSuHHjlJWVpZKSEi1durTJ8S+88IIWLlyovXv3auTIkRo8eLA2bdqkBx98UGvWrHHgE6A9mdhYmTET/G372fvOFQMAAIBOL6JnuIqKivToo4+qurpa99xzj84777xG+3fsaPwupi1btmjZsmVKTk7WI488ovT0dEm+0DZv3jwtWrRIw4cPV1JSUod9BrQ/M36i7KcrfY0tG2RLiqXUVGeLAgAAQKcU0TNcL7/8skpLSzVr1qwmYUuSBg8e3Kj91ltvSZJmzJjhD1uSlJWVpcmTJ6usrEwrV64MbtEIvqGjpC4pvu0aj+y6jxwtBwAAAJ1XxAaukpISffrpp0pISNBll1120uOrqqq0ZcsWSdKECROa7K/rW79+ffsWig5noqJkzmuweMYn7zlYDQAAADqziL2lcNu2baqurtbIkSPldru1evVqffXVV/J4POrTp4/OP/98paSk+I/Py8tTdXW1unTpoh49ejS53oABAyRJubm5HfUREETmgstkV/zN19iVLc/e3XKfkeloTQAAAOh8IjZw7d27V5LUtWtXPfDAA8rOzm60/+WXX9btt9+u888/X5JUWFgoSc2GLUmKi4tTYmKiSktLVV5ervj4+CBWj2AzZwyQzhgg7d0lSSpf+baSr7/d4aoAAADQ2URs4CotLZUkffDBB3K73br11ls1btw4VVRUaPny5Vq6dKl+85vfKCMjQ/3791dFRYUkKSYmJuA1Y2NjVVpaqoqKilMKXHPmzGnSFxMTo/nz50uSevbs2ZaPhnZSOvmbOv7sk5Kk8veXK+WGO5TK4hk4gdvt+zHK2EBzGB8IhLGBQBgb4Sdin+Hyer2SpJqaGt1www267LLL1KVLF/Xq1UvXX3+9JkyYII/Ho7/9zXdbmbVWkmSMCXjNumMQHuIuniJFRUmSvEUFqvwXy/4DAACgdSJ2hqtuBsoYo0suuaTJ/ksvvVSrV6/Wl19+2ej4ysrKgNesqqqS5Lu98FQsWLCgxf2FhYWEOKeNGCtt9AWtkhVvqToj09l6EHLqfgNZUFDgcCUIRYwPBMLYQCCMjdBjjGm0QnlrRewMV91gTklJUXR0dMD9x44dk1R/e9/hw4ebvV5FRYVKS0uVmJjI81thxHXBJP92xWerZMtKHawGAAAAnU3EBq66VQVLS0ubnUUqKSmRVD9blZGRoejoaBUXFzcbunbt8i2u0K9fv2CVDCeMGiclJfu2q6pk133obD0AAADoVCI2cPXr10+9evVSVVWVtm/f3mT/F198IUkaOHCgJN9iFiNGjJAkrV69usnxdX1jx44NVslwgHFHN34n18fvOlgNAAAAOpuIDVySdNVVV0mSnnvuORUXF/v7c3JytHTpUknS5MmT/f3Tpk2TJC1ZskT5+fn+/uzsbK1YsULx8fGn9BJldC7mwvrbCpWzTXY/71oDAADAqYnYRTMkadKkSdq8ebNWr16tu+++W1lZWaqsrNS2bdvk8Xg0adIkTZgwwX/8qFGjdMUVV+jtt9/Wvffeq5EjR6qmpkabNm2S1+vVnXfeqaSkJAc/EYLB9Bsk96Ch8uz8SpJkP/ynzLd/4HBVAAAA6AwiOnC5XC7dfffdeuedd/Tee+/5byMcNGiQJk+erIsvvrjJOTfeeKMyMzO1fPlybd68WVFRURoxYoSuueYaDR06tKM/AjpIwuQrVVwXuFa/L3vNDTLRgd/JBgAAAEiSsaw7HrLy8/NZFj5E9EiIV8G/Xylb6XsBtvnBPXKd1zSQI/KwfC9awvhAIIwNBMLYCD0sCw90AFdikuIaPMtlP3rHwWoAAADQWRC4gFMUP/nK+sbWjbKH8gMfDAAAAIjABZyy6KGjpPQz/G378QoHqwEAAEBnQOACTpExRuai+tcE2I/flfV4HKwIAAAAoY7ABbSCOf9SKap2cc9jRdKmNc4WBAAAgJBG4AJawSR3lRl7gb/tff/vDlYDAACAUEfgAlrJXHJ5fWPrRtkD+50rBgAAACGNwAW01pnDpT79/U27armDxQAAACCUEbiAVjLGNJrlsp+8K1tZ6WBFAAAACFUELqANzISJUmycr1FWIrvuQ0frAQAAQGgicAFtYOITfKGrlmXxDAAAADSDwAW0kZnYYPGM3dtld293rBYAAACEJgIX0Eam7wBp0FB/27631MFqAAAAEIoIXMBpMJdO82/btR/KFh9xsBoAAACEGgIXcBrM2AullO6+hscju+ofzhYEAACAkELgAk6DcbtlJl7hb9tVf5f1VDtXEAAAAEIKgQs4TebiqZI72tc4dkR23cfOFgQAAICQQeACTpNJ7ioz/mJ/2777lqy1DlYEAACAUEHgAtqBuezK+sbu7VLONueKAQAAQMggcAHtwPQbKGWd5W/bd99ysBoAAACECgIX0E5ck+pnuez6j2UPH3KwGgAAAIQCAhfQXs4eL/Xs7dv2emVXvOlsPQAAAHAcgQtoJ8YVJTP5Kn/bfvhP2dISBysCAACA0whcQDsyF35dSkz2NSorZD/gRcgAAACRjMAFtCMTGycz8XJ/2777Fi9CBgAAiGAELqCdmcumSW63r3GsSPazD5wtCAAAAI4hcAHtzHTpJnP+Zf62/efrvAgZAAAgQhG4gCBouHiG8vZIm9Y5VwwAAAAcQ+ACgsCknyGNPs/f9v79r8xyAQAARCACFxAkrsv/rb6x8ysp+wvnigEAAIAjCFxAkJhBQ6UhI/1t79t/dbAaAAAAOIHABQSR64pr6xtffi67e7tzxQAAAKDDEbiAYBo2Wso809/0/v1VB4sBAABARyNwAUFkjGk8y7XhU9m8Pc4VBAAAgA5F4AKCbfR5UvoZ/qblWS4AAICIQeACgsy4XDINZrnsmg9lD+xzsCIAAAB0FAIX0AHMeV+TevfxNaxXdtkrzhYEAACADkHgAjqAcUXJTJ/pb9vPPmCWCwAAIAIQuIAOwiwXAABA5CFwAR2EWS4AAIDIQ+ACOpA572tSWoNZrqV/cbYgAAAABBWBC+hAvlmub/vbds0Hsvt5LxcAAEC4InABHcyce1H9e7mslfeNF50tCAAAAEFD4AI6mHFFyXX19+o7/rVadle2cwUBAAAgaAhcgBPOOV/qP9jf9L7+RweLAQAAQLAQuAAHGGPkmjG7vmPrRtmtG50rCAAAAEFB4AKcMuxsachIf9P7+h9lrXWuHgAAALQ7t9MFhJKSkhLdfffdKi4uVkZGhp544omAx65atUrLly/Xvn375Ha7lZWVpRkzZmjIkCEdVzA6NWOMXN+aLe/8e30du7Klzz+VxlzgbGEAAABoN8xwNfDCCy/o+PHjp3TcwoULtXfvXo0cOVKDBw/Wpk2b9OCDD2rNmjUdUCnChRk0VBp9nr/tfe0Psh6PgxUBAACgPRG4am3evFmrVq3SpEmTWjxuy5YtWrZsmZKTk/XLX/5S9957r+bOnauHHnpILpdLixYtUklJSQdVjXDgmnG9ZGr/Kh7Kk/1gubMFAQAAoN0QuCRVVVXpmWeeUd++fXXllVe2eOxbb70lSZoxY4bS09P9/VlZWZo8ebLKysq0cuXKoNaL8GIy+sl8bbK/bd/6s2xZqYMVAQAAoL0QuCT99a9/1cGDB3XzzTcrKioq4HFVVVXasmWLJGnChAlN9tf1rV+/PjiFImyZb35Xio3zNUqKZZe/6mxBAAAAaBcRH7hyc3O1dOlSTZw4UcOHD2/x2Ly8PFVXV6tLly7q0aNHk/0DBgzwXxNoDdO1m8zUGf62XfGWbFGBgxUBAACgPUR04PJ6vXr66aeVkJCgWbNmnfT4wsJCSWo2bElSXFycEhMTVVpaqvLy8natFeHPTLla6trd16iukuVlyAAAAJ1eRC8Lv3z5cu3YsUO33367kpOTT3p8RUWFJCkmJibgMbGxsSotLVVFRYXi4+NbvN6cOXOa9MXExGj+/PmSpJ49e560JnQMt9v3VyU1NTWo36ds1v9T8cJHJUl29fvq+q3vKWbIiKB+T5yejhob6JwYHwiEsYFAGBvhJ2JnuAoLC/XnP/9Zw4cP18SJE0/pnLqX0hpjTnoM0Bbxl02TO/NMf/v44gWyXq+DFQEAAOB0ROwM1+LFi+XxeHTzzTef8jl1M1aVlZUBj6mqqpLku73wZBYsWNDi/sLCQgJciKj7LVNBQfCfq7LXfl/65U8kSdXbt+rQW6/IdUHLryuAczpybKDzYXwgEMYGAmFshB5jTKPVyVsrYgPXhg0blJiYqMWLFzfqr66uluQLO/PmzZMk3XfffYqLi/Pf4nf48OFmr1lRUaHS0lIlJiae9HZCIBCTNUJm3EWy6z6SJNklf5Adc75MXILDlQEAAKC1IjZwSVJpaam+/PLLZvdVVVX599XU1EiSMjIyFB0dreLiYh0+fLjJ4hm7du2SJPXr1y+IVSMSmH+7UXbjGqm6Sjp2RHbZX2WuucHpsgAAANBKERu4XnnllWb7Dx06pB/+8IfKyMjQE0880WhfTEyMRowYoc8//1yrV6/WtGnTGu1fvXq1JGns2LFBqRmRw/ToJTN1huzSP0uS7Dt/k71wkkxaX4crAwAAQGtE7KIZbVUXspYsWaL8/Hx/f3Z2tlasWKH4+HhddtllTpWHMGK+cY3UvXalyhqPvC8/zTN9AAAAnUzEznC11ahRo3TFFVfo7bff1r333quRI0eqpqZGmzZtktfr1Z133qmkpCSny0QYMLGxcl13s7y/9b0mQFs3yq77WObci5wtDAAAAKeMwNUGN954ozIzM7V8+XJt3rxZUVFRGjFihK655hoNHTrU6fIQTsacL40YI23ZIEmyryyWHTmGBTQAAAA6CWO5Rylk5efncwtZiHByiVZ7KE/eB++UPL4VNM3kq+S67qYOrwPNY/letITxgUAYGwiEsRF6TndZeJ7hAkKc6ZUhc/k1/rZ99y3ZfbscrAgAAACnisAFdALmG9dIqWm+htcr7x8WynprnC0KAAAAJ0XgAjoBExMr13dvre/YlS373jLnCgIAAMApIXABnYQZMUZm/CX+tn3jRdnDhxysCAAAACdD4AI6ETPzZikp2deorJD3pd+ysAoAAEAII3ABnYhJ7ipz3c31HZvXya75wLmCAAAA0CICF9DJmAkTpeFn+9v2z8/IFh91qhwAAAC0gMAFdDLGGLlm3S7FxPo6Soq5tRAAACBEEbiATsikpslcc0N9x4ZPZNd+6FxBAAAAaBaBC+ikzMQrpCEj/W378tOyx444VxAAAACaIHABnZRxueS64U4pNs7XUXpc3hcXcWshAABACCFwAZ2YSU2T+bcb6zv+9Znsp+85Vg8AAAAaI3ABnZy5+BvSsNH+tn35d7IFBxysCAAAAHUIXEAnZ1wuuW78Dykh0ddRWS7vs4/L1tQ4WxgAAAAIXEA4MN1TZWbdXt+xY6vs8tecKwgAAACSCFxA2HCd+zXfS5Fr2bf+JLtru3MFAQAAgMAFhBPznf8ndU/1NWpq5H3ml7LlZc4WBQAAEMEIXEAYMQmJct00RzK1f7ULDsj+cSFLxQMAADiEwAWEGZN1lsz0mf62Xfuh7EfvOFgRAABA5CJwAWHITL9OGjLS37Z/+p3s/lwHKwIAAIhMBC4gDBlXlFw3z5GSuvg6qqvkffoXshXlzhYGAAAQYQhcQJgyKT3k+ve76zvy98r+4f94ngsAAKADEbiAMGZGjpP5xjX+tl37oex7yxysCAAAILIQuIAwZ66e1fh5rr/+XnbnVw5WBAAAEDkIXECYM1FRct1yj5TS3ddRUyPvb38uW3zE2cIAAAAiAIELiACmSze5/t+9UlSUr+PoYXmf+rmsp9rZwgAAAMIcgQuIEGbwcJlr/72+Y8eXsn9+xrmCAAAAIgCBC4gg5rLpMhdM8rftquXyvv93BysCAAAIbwQuIIIYY2Rm3SYNyPL32T//TjZ7i4NVAQAAhC8CFxBhTHSMXLf/WOraYBGNpx6VPZTnbGEAAABhiMAFRCCT0sMXutzRvo6S4/L++mHZ0hJnCwMAAAgzBC4gQpmBQ2Ru/I/6joP7fTNdrFwIAADQbghcQARzjb9E5qrv1nds2yz7x0Wy1jpXFAAAQBghcAERzkybKTNhor9tP3lX9s0/OVcQAABAGCFwARHOGCNz/Z3SmcP9fXbpn+VdtdzBqgAAAMIDgQuATHS0XHfMldLP8PfZl34r+/lqB6sCAADo/AhcACRJJjFZrrvmSSk9fB3WK+8zj8lu/9LRugAAADozAhcAP9MjVa67HpTiE30d1VXy/uZ/ZHN3OlsYAABAJ0XgAtCI6Zvpu72w7h1d5WXyPvGgbP4+ZwsDAADohAhcAJowQ0bIdet9UlSUr6OkWN4FP5UtPOhsYQAAAJ0MgQtAs8zoc2W+f7dkjK/j6GFf6CoqdLQuAACAzoTABSAg1/hLZL53W31HwQF5f3W/7NHDzhUFAADQiRC4ALTIdck3ZK79fn3HoTxf6Dp2xLmiAAAAOgkCF4CTck35lsyMG+o7Duyvnekqcq4oAACAToDABeCUuC6/RubqWfUd+Xvl/cV9LKQBAADQAgIXgFPmmnadzDe/W99RcEDeX/xY9gBLxgMAADTH7XQBTqmsrNTGjRu1fv167dy5UwUFBfJ6vUpLS9P48eM1ffp0xcXFNXvuqlWrtHz5cu3bt09ut1tZWVmaMWOGhgwZ0sGfAuh4riu/LW9MjOyrz/s6jhTK+4sfy/Wf/yNzxgBHawMAAAg1ETvD9dFHH+mxxx7TypUrZa3V6NGjNXToUB06dEivvPKKfvzjH+vYsWNNznvhhRe0cOFC7d27VyNHjtTgwYO1adMmPfjgg1qzZo0DnwToeK6pM2Rm3V6/ZPzxY/I+9hPZnV85WxgAAECIidgZLrfbrSlTpmjatGlKT0/39x85ckTz58/Xrl279Pzzz+uuu+7y79uyZYuWLVum5ORkPfLII/7zsrOzNW/ePC1atEjDhw9XUlJSh38eoKO5LvmGvLFxss89IXm9UlmpvI8/INcdc2WGjXa6PAAAgJAQsTNcl1xyiW6++eZGYUuSunXrpptuukmStGbNGnk8Hv++t956S5I0Y8aMRudlZWVp8uTJKisr08qVKzugeiA0uCZMlOvW+yR37e9uKivk/fX/yG5kthcAAECK4MDVkv79+0uSqqurdfz4cUlSVVWVtmzZIkmaMGFCk3Pq+tavX99BVQKhwZwzQa47fyrFxPo6PNXyLvyZvO+/7WxhAAAAIYDA1YyDB33LXEdFRflvD8zLy1N1dbW6dOmiHj16NDlnwADfYgG5ubkdVygQIszwc+T6z4ek+ARfh/XKvvRbeV/5vay3xtniAAAAHETgasbbb/t+M3/22WcrOjpaklRYWChJzYYtSYqLi1NiYqJKS0tVXl7eMYUCIcQMHi7Xfz0qpdT/HbHv/E3ep34uW1nhYGUAAADOidhFMwLZsGGDVq5cqaioKM2cOdPfX1Hh+wdjTExMwHNjY2NVWlqqiooKxcfHn/R7zZkzp0lfTEyM5s+fL0nq2bNna8tHkLhrn1FKTU11uJIQl5qqml89qyP/+1/y5GT7+v61Wq4nHlS3ub9UVLfmf2HRmTE20BLGBwJhbCAQxkb4YYargX379uk3v/mNrLWaPXu2MjMz/fustZIkU7cMdjPqjgEiWVSPXur+v08pdtwF/j7Pjq06fO/Nqs7d6WBlAAAAHY8ZrlqHDx/Wz372M5WWlmr69Om64oorGu2vm7GqrKwMeI2qqipJCvjC5BMtWLCgxf2FhYWEuBBR91umgoIChyvpPOwP/ksmebHsymWSJG/BAR2+7xa5bvkvmRFjHa6u/TA20BLGBwJhbCAQxkboMcY0Wdm8NZjhklRcXKxHHnlEhYWFmjhxombPnt3kmLrb+w4fPtzsNSoqKlRaWqrExMRTup0QCHfGFSXXd/+fzLd/UP+C5PIyeX/9P/Iu/Yus1+tsgQAAAB0g4gNXeXm5Hn30Ue3fv1/nnXeebr311mZvG8zIyFB0dLSKi4ubDV27du2SJPXr1y/oNQOdiWvSlXLd/pP6ZeOtlf3bS/Iu+plsWYmzxQEAAARZRAeu6upq/eIXv9DOnTs1evRo3X333XK5mv+fJCYmRiNGjJAkrV69usn+ur6xY8PnVimgvZizx8v1419IqWn1nRvXyPu/P5Ldt9uxugAAAIItYgOX1+vVk08+qS+++ELDhg3TPffc418VJpBp06ZJkpYsWaL8/Hx/f3Z2tlasWKH4+HhddtllQa0b6KxM3wFy3b9AGnVufeehfHkfvUfez1Y5VxgAAEAQReyiGcuXL9eaNWskScnJyVq8eHGzx82ePVtdunSRJI0aNUpXXHGF3n77bd17770aOXKkampqtGnTJnm9Xt15553+FyUDaMokJMl1x1zZt/8q++bLkrVSVZXs4l/Jm7NN5trvy7ijnS4TAACg3URs4CopqX92pC54Nefaa6/1By5JuvHGG5WZmanly5dr8+bNioqK0ogRI3TNNddo6NChQa0ZCAfG5ZKZPlM2c7C8z/xKqn2Oy763VHbHVrl+8COZtL4OVwkAANA+jGXd8ZCVn5/PsvAhgiVag8MWHJD3t/OlPTn1nTGxMjNvlvnalBbfexcqGBtoCeMDgTA2EAhjI/SwLDyATsukpsn13z+Xufgb9Z1VlbJ/XCjvU4/KlhQ7VxwAAEA7IHABcJSJiZVr9u1y3fETKSm5fsfnq+V96D9kv/yXY7UBAACcLgIXgJBgzp4g14O/loafXd95tEjexx+Q96/PylZVOlYbAABAWxG4AIQMk9JDrrvmyVz771KD1zTYf74h70N3yWZ/4WB1AAAArUfgAhBSjMsl15Sr5frxY1L6GfU7DuXJ+8sfy/vSb2UrypwrEAAAoBUIXABCkuk3UK65C2SmXC2Z+h9V9v235X3wh7Kb1ztXHAAAwCkicAEIWSY2Vq5r/12uH/9C6tO/fkdRoby/fkje3z/OSoYAACCkEbgAhDwzIEuu+xfIXPkdKarBs12rV8r709vk/eAfst4aBysEAABoHoELQKdg3NFyffM7cv30cSnzzPodJcd97+362X/J7vzKuQIBAACaQeAC0KmYPv3l+vEvfCsZxsbV78jdIe/8e+V9/knZ4iPOFQgAANAAgQtAp2NcUb6VDB9+Sua8Sxrtsx+/K+/9t8m74k1Zj8ehCgEAAHwIXAA6LdOth1w/+JFc//UzqW9m/Y7yMtm/LJZ33p2yGz6RtdaxGgEAQGQjcAHo9EzWCLnuf1zmO7dICYn1Ow7ul/ep+fI++l+y2VucKxAAAEQsAheAsGCiouS6bLpcj/xW5mtTGr27S7uy5f3lT1Tzm4dl9+c6VyQAAIg4BC4AYcUkd5Xr+h/K9eCvpdHnNd65aa28D/2HvM89KVtwwJkCAQBARHGf/BAA6HxMn36K+uH9stu/lPe156W6JeOtlf3kXdnVK2XOv1TmimtlemU4WisAAAhfzHABCGvmzOFy/ffP5br9J1Ja3/odXm/tioa3y/v7x2UP7HeuSAAAELaY4QIQ9owx0jkT5Bp1rm926+2/SoUHfTutV3b1StnPVsmc+zWZadfKZPRztmAAABA2CFwAIoaJipL52hTZ8y+T/WyV7LK/SHXPclmv7JpVsmtWSaPOlWvKt6Sss3xhDQAAoI0IXAAijnG7ZS6cJDthouyaD2SXvSIdbHBL4aa18m5aK/UfLDPlapmxF8pERTlXMAAA6LR4hgtAxDJRUXKdf6lc//N/Mjf/SDrxVsLcHbLPPCbvT26R952/yZaXOVMoAADotJjhAhDxjCtKZvwlsuddLH2xQd5/viFt3Vh/QFGB7Cu/l/3byzITLpGZeLlM3wFOlQsAADoRAhcA1DLGSCPGKmrEWNk9ObLv/E127QdSTY3vgMpy2VXLZVctlwYPk5l4hezUb8pExzhaNwAACF3GWmudLgLNy8/PF/95QkNqaqokqaCgwOFK0NFsUaHse0tlP/yHVFbaZL+razfFf326ysdcJNMr3YEKEcr42YFAGBsIhLEReowxSk9v+//HE7hCGIErdPDDD7ayUnbdh7Ir35ZydzR/UNZZMhdM8i2yERffsQUiJPGzA4EwNhAIYyP0nG7g4pZCADgFJjZW5sKvSxd+XXbXdtn335Zd+6FUXVV/UPYXstlfyP7pd77QdeEk6UyWlgcAIJIRuACglcyAM2UG3CV73b8rceNnKnvnTdXs211/QGWF7wXLn7wrpabJnHex74sXKgMAEHG4pTCEcUth6GB6H4GkpqbKWquCNR/7QtaaD6Xyps96SZL69PcFr3O/JpOa1rGFwhH87EAgjA0EwtgIPdxSCAAOM8bIDBwiM3CI7HU3yX6+WvaT96St/5Ia/tJkf67s63+Uff2P0oAsmfO+JnPO+TI9ejlWOwAACC4CFwC0IxMTKzP+Emn8Jb4VDtd96Jv1OnGhjV3ZsruyZf/ye6nfIJlzJsiMOV9KP4NnvgAACCMELgAIEtO9p8yUb0lTviV7YH99+Mrf2/jAPTtl9+yU/dtLUu8+vvB19nhpwJkyrihnigcAAO2CwAUAHcCk9ZGZ/m3ZaTOl/btl13wo+/mn0oH9jQ88uF92+Wuyy1+TkrrIjBgjjRgrM2KMTGKyM8UDAIA2I3ABQAcyxkh9B8j0HSDNuF42f6/shk9lP1/d9LbDkmLZ1e9Lq9+XNS5p0FCZUeNkhp8jnTFAxuVy5DMAAIBTR+ACAAeZ9DNkpp0hTbtO9vAh34Ib//pM2vGlVFNTf6D1Sju+lN3xpeySP0iJydLQkTJDR8sMHy2lpvPsFwAAIYjABQAhwvToJfP1b0pf/6ZsWam09V+ym9bJbl4nHT/W+ODS49L6T2TXfyIrSd1TZYaNloaNlhk6SqZrNyc+AgAAOAGBCwBCkElIlMZeKDP2Qlmv17ewxqZ1sl9skHZvl7zexicUFch+vEL6eIUvgPXpL5N1ljR4uMzg4TLdezrxMQAAiHgELgAIccblkjLPlMk8U/rmd3yzX9u/kN26UXbrRilvT9OT9ufK7s+VVr7tC2A9eskMGiadOUxm8HApox/PgAEA0AEIXADQyZiERGn0eTKjz5Mk2aNFsl9tkrZulP1qo1RU2PSkw4dkDx+S1qzyBbD4RN8iHIOHyQwaKvUfLBOf0KGfAwCASEDgAoBOzqR0l5kwUZowUdZa6VC+L4Dt2Cq740up8GDTk8pLpS3rZbes9wUwY3zvAMs8s3Y2bbDUb6BMdEwHfxoAAMILgQsAwogxRuqdIdM7Q7rkG5Ike/RwbfjaKrv9S2nvLt+qhw1ZKx3YJ3tgn7R6pS+ERUX5ngWrC2H9B0np/WSiozv8cwEA0FkRuAAgzJmUHtK4i2TGXSRJshVlUk62b4n5HVt97/8qK216Yk2NtCdHdk+O9ME/6kNYWl+ZMwb43id2xgDfO8GSu3boZwIAoLMgcAFAhDFxCdLws2WGny1JvlUQCw7I7sqWdm+X3b1d2psjVVU1Pbmmpn5BDr3vC2GSlNK9NoBlSn0yZfr0k3r3ZTYMABDxCFwAEOGMy1V/G+KEiZIkW1Mj5e3xha+6ELZ/j1Tjaf4iR4uko0WyW9b7zvddWOqVLmWcIZPez/dnRj8prQ/PhgEAIgaBCwDQhImK8t0qeMYA6WtTJEnW4/E957V3l7Rvl+/PvbukkuLmL2K90sH90sH9sp+v9nVJviCWmuYLYL0zpF4ZMr37SL0zpK7dfM+hAQAQJghcAIBTYtxuqW+mTN9MSZdKkm9VxGNF0t7dsvt8Aczm7ZEO7A88G2a90qE86VCe/5ZE/62JsfFS73SZXhm+ANarduYtrY9MYnJwPyAAAEFA4AIAtJkxRkrpIaX0kBk51t9vPR6p4IDvtsT8PVLeXl8QO7hf8gQIYpJUWV6/UEfdteo2EpN9tz72TJN69pZ69pLp2du33T3VNysHAECIIXABANqdcbul9L5Sel8ZXeDvtzU1UkG+L4gd2C8dzJM9lCcdzJOOH2v5oqXHpZxtsjnb6q9Xt+FySd16Sj1714ewhtvcqggAcAiBqw2qqqr0xhtv6OOPP1ZhYaGSkpI0evRozZw5Uz169HC6PAAIWaZ2WXml9dWJ8ceWlUgH82UP7vfdcngwT/ag79ZDlZe1fGGvVzp8SDp8SHbb5vpr1m24o6XuPaVuPWW695S6pUrde8p0T63vT0hsz48KAIAkAlerVVVV6eGHH9a2bdvUrVs3jRs3TgUFBXr//fe1YcMGPfLII0pLS3O6TADodExCkjTgTJkBZzbqt9ZKx49KB2pnwwoPSoUHZQsPSoWHfM+QnYynWjqULx3Krw9hUqNtxcX7ZskaBrGUHjIp3X3L3nftLiUm+1Z1BADgFBG4Wun111/Xtm3blJWVpfvvv19xcXGSpKVLl+oPf/iDnnrqKT300EMOVwkA4cMYI3XpJnXpJpN1VpP9tqpSOlxQH8IONwhjhQd9tyKeiopyKX+vlL83cCiLcktdu/m+Urr7wlhXXyAzXRsEs6RkbmEEAEgicLWKx+PR8uXLJUk33XSTP2xJ0vTp07Vq1Spt3bpVOTk5GjhwoFNlAkBEMTGxDZ4Xa8qWl/luNzxSKHu4QDpSKBUVyh4plIpq2y0t5NFQjcd3TlGB79oNv0/D49xuqWt3He7ZS66U7vLGxktdUqTkrlJyikyXFKlLV19fQhLhDADCGIGrFb766iuVlpaqd+/eGjBgQJP948ePV25urtatW0fgAoAQYeITpL6ZviXtm9lvvV6p5JhU1DCIFfoC2tHDvpc6HyuSqqpO/Zt6PNLhQ6o+fKjZ3Y1nzaJqg5gvgJnklPowltzVF86SU2rbXWTc0adeBwDAcQSuVsjNzZWkZsOWJH/IqjsOABD6jMvlv2VRmWc2H8qs9S3ccaxIOlokWxfCjta2jxVJx4742tWtCGaSVFPjv450Qhhrpq34RCkpWUrqIiV1kUlKlhK7+PtMbb//mMRk36qRAABH8BO4FQoLCyUp4EqE3bt3b3QcACA8GGOkhETfV/oZzYYyqTaYlZX6Z8WSvdXyHi1SSX6eVHxU9vhRqfiobwn848d8Yau1ykt9XwUHfN/zxBqaOyc+wfces4YhrTaMNWonJUsJyVJikhQTy62OANAOCFytUFFRIUmKjY1tdn/dM111x53MnDlzmvTFxMRo/vz5kqSePXu2pUwEgbv2t8OpqakOV4JQw9hA83x3QtSNj8RmnhGzXq9s6XF5jx2R92iRvMeOqKb2T++Jfx474nsWra3Ky3xfhQd93/vEWpo7xx0tV1KyTFKyXInJciV3kUnsUt+XVLud6Ntu1BfT/P9Poh4/OxAIYyP8ELhawdpm/y/plPcDAFDHuFwyyV3lSu7qe8bsJGxlhWqOHpE9flTe48fkLT4mb3H9tj1e1y6W9/hReYuPtf72xoY8vtk5HS1Sq+fhYmLkSuxSG9KSmwlpdX8myZWYLJOQKFdikkxiskxcPEvvAwgrBK5WiI+PlyRVVlY2u7+uv+HqhS1ZsGBBi/sLCwsJcSGi7rdMBQUFDleCUMPYQEvafXy43FLXnr6vFhhJLmulqkqp5LhUUiyVFsseL/a1S4t9fSXHZUuKfUvnl9T2tWZxkECqquStKpT3SBtusTdGikvw3b4Zn+D7SkjyLX4SnyDFJ0kJ9dumwTG+PxOl6JiQvx2Snx0IhLEReowxSk9Pb/P5BK5WqLvF7/Dhw83uLyoqanQcAABOMcZIsXG+rx6+f8CdSgSx1dVSWYnvq9T3ZRu2ywL3yVN9+oVbW/+cWsPuQIc31xnlrg9f8Yn+bV84q+1L8G2buv1xCVJ8fO2fCTzDBqDdELhaoX///pKkXbt2Nbs/Jyen0XEAAHQ2Jjq6/uXOdX2neK6tqqwNX6W+WbOyE4JZbTizZbX7y8t8waqs1Dcb115qPPUzdg3ra67mQNcwrvoAFhdfG8riZeoC2YkBLS5BplG7/lwTFdV+nw1Ap0PgaoWhQ4cqISFBBw8e1K5du5osD//ZZ59JksaMGeNEeQAAOMrExEoxsVJK/Wq+pxzWPB6poswXvhoEMVteJpWXSGW1C3+Ul9T21R1bWn9OzSm+wPqUCvL6rlt2ajNtLe6Lia0PabXh7UjXFJmERHnlOiHAtRDq3NHMugGdEIGrFdxut77xjW9oyZIlevbZZzV37lz/81pLly5Vbm6uhg4dqsGDBztcKQAAnYtxu/3L1jfqP8XzrbW+RUICBbbystrQ5tu25WW+mbe68FZR7gt8bVmq/2SqKn1fx474u1qazwsY3KLctTNnDb5i60Laif1xjftjG8zW1X4x8wZ0DAJXK82YMUObN2/Wtm3bdNddd2no0KEqLCzU9u3blZycrNtvv93pEgEAiDjGGN9MUkyslNK9vr8V17DW+p5DKy/zha/y2hBWXiZ7Qrtuv61o2Ndguz1vkaxT4/Hdill6vGntgT5TS9eLjmkUzvxBLDa+abCrDW1Ng13ddpyMiwAHNIfA1UoxMTF68MEH9frrr+ujjz7S2rVrlZiYqEsuuUQzZ85kwQwAADopY4wvhETHSF1SGu9r5bVsTU3jAHZCQEuKcslbVqqyw4WNQ11F+QmBr9x3e2MwVFf5vo4fa1x7S5+rpevFxDYT4BJOmGWLPyHcxTUf7mLieD0AwgaBqw1iYmI0c+ZMzZw50+lSAABACDJRUVJiku+rYX/tn4m1S39XnGTpb1u3vP+Js26VdeGt+S9bWVHfrmyw73TezXYydbdOnvgZWvp8LV2vYUBrEMxMw3asb3bNfwtlbDP7ao837ujT/YRAmxC4AAAAQlSj5f3VvfG+NlzPejxSgDBmTwxnLfXXtT3tuFDJiSprv0/jCbi2B7god4MwFtcozJnYuPqAFxvXIKg1mIU7IdwplufgcGoIXAAAABHCuN2Su+nMm9TWAFfddIatNtDVz8CVNQ12DY5rFOKCsWhJnRpP/TvmTvwcLX3Glq7pjm42wPlCXPMBrtG+EwKcYmNP91MiBBG4AAAA0CbGHS0lRTdZXVJqw3NvdYuWNDeTdmJQ8/dX1N8+WdkgxNUdUxXEWyglX70l1U3e+Sa1PcQdiImVKz5B3pjYxjNysQ1vpzxhRq6lfTGxPA/nMAIXAAAAHNdo0ZLkrk33t+GatqamNnw1CGmVFb6gVlFWv69RYKuov43yxAAX7NsoJamqUt4Aq1y2KcQZI8XEnTCbFtf4dslGz8WdPNwpJoZ3wrUCgQsAAABhyURFSQmJvq8T97XxmtZTLVVWnhDgymsXMmkuwNXNxAUIcBXlkjdIK1FKkrX1z8OduKul01q6pnHVB7jYE26lPOHWyka3Uvpvo2y4P8HXF8Yv9iZwAQAAAKfIuKN9z26113Nw1vpmzWpDWreEeNmKch09kNfyLZMNV6NsZgYvaK8TkHzXrntp+Im7WjqtpWu6XPVB7IQZNxMbJ3P+pTIjx51u5Y4gcAEAAAAO8d1KGe37Su6i6NpXBpjuvX1/tuGavtcJVJ0Q0nwzbaoMcMtkw1cKNNlX0ewMWbvyeqXyUt/XiZ9HkgYNbfOspNMIXAAAAEAY8b1OIDbgqodtCnFer+89a6cU4Gr3BQxwtX8GeFatWbFxbag6NBC4AAAAALTIuFz1C2x07dZ4Xxuvab01vtBVccKtkZUVtc/D1Qc0kzn49D+EQwhcAAAAADqccUXVLpqR0HSfA/UEC4vyAwAAAECQELgAAAAAIEgIXAAAAAAQJAQuAAAAAAgSAhcAAAAABAmBCwAAAACChMAFAAAAAEFC4AIAAACAICFwAQAAAECQELgAAAAAIEgIXAAAAAAQJAQuAAAAAAgSAhcAAAAABAmBCwAAAACChMAFAAAAAEFC4AIAAACAICFwAQAAAECQELgAAAAAIEgIXAAAAAAQJAQuAAAAAAgSt9MFIDBjjNMl4AT8N0EgjA20hPGBQBgbCISxETpO97+FsdbadqoFAAAAANAAtxQCAAAAQJAQuIBTcN999+m+++5zugyEIMYGWsL4QCCMDQTC2Ag/PMMFnIKqqiqnS0CIYmygJYwPBMLYQCCMjfDDDBcAAAAABAmBCwAAAACChMAFAAAAAEFC4AIAAACAIOE9XAAAAAAQJMxwAQAAAECQELgAAAAAIEgIXAAAAAAQJAQuAAAAAAgSAhcAAAAABAmBCwAAAACChMAFAAAAAEHidroAIFy8+uqreuWVVyRJd911ly688EKHK4IT9u/fr7Vr12rTpk3Kz8/XsWPHlJiYqCFDhmjatGkaNmyY0yUiyKqqqvTGG2/o448/VmFhoZKSkjR69GjNnDlTPXr0cLo8OKSyslIbN27U+vXrtXPnThUUFMjr9SotLU3jx4/X9OnTFRcX53SZCBElJSW6++67VVxcrIyMDD3xxBNOl4TTQOAC2kFeXp5ef/11GWPEu8Qj28MPP6yioiLFx8frzDPPVFZWlvbt26c1a9Zo7dq1uv766zVt2jSny0SQVFVV6eGHH9a2bdvUrVs3jRs3TgUFBXr//fe1YcMGPfLII0pLS3O6TDjgo48+0tNPPy1JOuOMMzR69GiVl5crOztbr7zyij7++GPNmzdPXbt2dbhShIIXXnhBx48fd7oMtBMCF3CarLV6+umnlZCQoDPPPFPr1q1zuiQ4qG/fvpo1a5YmTJggt7v+R+w777yjZ555Rn/84x81evRo9e3b18EqESyvv/66tm3bpqysLN1///3+GYulS5fqD3/4g5566ik99NBDDlcJJ7jdbk2ZMkXTpk1Tenq6v//IkSOaP3++du3apeeff1533XWXg1UiFGzevFmrVq3S17/+da1YscLpctAOeIYLOE3vvvuutm7dquuvv16JiYlOlwOH3X///brooosahS1Jmjx5skaPHi2v16tPP/3UoeoQTB6PR8uXL5ck3XTTTY1uD5s+fbr69++vrVu3Kicnx6kS4aBLLrlEN998c6OwJUndunXTTTfdJElas2aNPB6PE+UhRFRVVemZZ55R3759deWVVzpdDtoJgQs4DUePHtVLL72kkSNH6mtf+5rT5SDE9e/fX5LvN9oIP1999ZVKS0vVu3dvDRgwoMn+8ePHSxKz4Gii7mdDdXU1t5FFuL/+9a86ePCgbr75ZkVFRTldDtoJgQs4Dc8++6yqqqp08803O10KOoGDBw9KklJSUpwtBEGRm5srSc2GLUkaOHBgo+OAOnU/G6KiopSUlORwNXBKbm6uli5dqokTJ2r48OFOl4N2ROAC2mj9+vVavXq1vvWtbzW5RQQ40YEDB7RhwwZJ0rhx4xyuBsFQWFgoSQFXIuzevXuj44A6b7/9tiTp7LPPVnR0tMPVwAler9f/PPisWbOcLgftjMAFtEFFRYUWL16s9PR0XXXVVU6XgxBXU1OjRYsWqbq6WhdccIF/pgPhpaKiQpIUGxvb7P66Z7rqjgMkacOGDVq5cqWioqI0c+ZMp8uBQ5YvX64dO3Zo9uzZSk5OdroctDNWKURE+tWvfqW9e/e26pwf/vCHGjx4sCTp5Zdf1uHDh/XAAw/w28gwc7pjoznPPvusvvrqK/Xu3ZvbT8PYyV4JwSsjcKJ9+/bpN7/5jay1mj17tjIzM50uCQ4oLCzUn//8Zw0fPlwTJ050uhwEAYELEamgoEB5eXmtOqeyslKStGPHDv3jH//QxRdfrBEjRgSjPDjodMZGc1599VW988476tq1q+bOncvzGWEsPj5eUuDxUNfPy20hSYcPH9bPfvYzlZaWavr06briiiucLgkOWbx4sTweD7+QC2MELkSk+fPnt/ncDRs2yFqrPXv2aN68eY327d+/X1L9P7InTJigb3zjG6dTKjrY6YyNEy1fvlyvvPKKEhISNHfuXF54G+Z69uwpyfcP6eYUFRU1Og6Rq7i4WI888ogKCws1ceJEzZ492+mS4KANGzYoMTFRixcvbtRfXV0tyTcDVvfvjfvuu49f2nRCBC6gjXbv3h1w3/79+7V//35uD4lgH374oZ577jnFxsbqvvvuYyxEgLqlvXft2tXs/rr3b9Udh8hUXl6uRx99VPv379d5552nW2+9VcYYp8uCw0pLS/Xll182u6+qqsq/r6ampiPLQjshcAGtdN111+m6665rdt/ChQu1atUq3XXXXbrwwgs7uDKEig0bNmjRokWKiorSPffco6FDhzpdEjrA0KFDlZCQoIMHD2rXrl1Nlof/7LPPJEljxoxxojyEgOrqav3iF7/Qzp07NXr0aN19991yuVi/LNK98sorzfYfOnRIP/zhD5WRkaEnnniiY4tCu+JvOQC0o6+++koLFiyQJN19990aPXq0wxWho7jdbv8txM8++2yj1QiXLl2q3NxcDR06tMUFVhC+vF6vnnzySX3xxRcaNmyY7rnnHrnd/N4biAT8TQeAdvTzn/9cVVVV6tWrl9auXau1a9c2OWbo0KGaNGmSA9Uh2GbMmKHNmzdr27ZtuuuuuzR06FAVFhZq+/btSk5O1u233+50iXDI8uXLtWbNGklScnJyk+d16syePVtdunTpyNIABBmBCwDaUWlpqSTfrSCHDh0KeByBKzzFxMTowQcf1Ouvv66PPvpIa9euVWJioi655BLNnDmTBTMiWElJiX+7Lng159prryVwAWHGWF4MAgAAAABBwTNcAAAAABAkBC4AAAAACBICFwAAAAAECYELAAAAAIKEwAUAAAAAQULgAgAAAIAgIXABAAAAQJAQuAAAAAAgSAhcAAAAABAkBC4AAAAACBICFwAAAAAECYELAAAAAIKEwAUAAAAAQULgAgAAAIAgIXABAAAAQJAQuAAAAAAgSAhcAAC0gzfeeEPXXXedvvOd72jHjh3NHrNhwwbNnDlT1113nT766KMOrhAA4AQCFwAA7eCqq67SyJEjVVNToyeffFLl5eWN9h85ckSLFi2StVYXX3yxLrroIocqBQB0JAIXAADtwBijO++8U127dtXBgwf1zDPP+PdZa/V///d/Ki4uVlpamm6++WYHKwUAdCQCFwAA7SQlJUW33367jDH66KOP9P7770uS/va3v2nz5s2KiorSXXfdpbi4OGcLBQB0GAIXAADt6JxzztG0adMkSc8++6w++OAD/eUvf5Ekfec739GgQYOcLA8A0MGMtdY6XQQAAOHE4/Ho/vvvV05Ojr9v9OjR+slPfiJjjIOVAQA6GjNcAAC0M7fbrdtvv93fTkhI0B133EHYAoAIROACACAIVqxY4d8uLy/X7t27nSsGAOAYAhcAAO1s/fr1Wr58uSSpf//+stZq4cKFOnr0qLOFAQA6HIELAIB2VPe+LUmaOHGiHnroIaWmpurYsWNauHCheHQaACILgQsAgHbi9Xr1f//3fzp+/LjS09P17//+70pISNBdd92lqKgobdy4UUuXLnW6TABAByJwAQDQTt58881m37eVlZWlf/u3f5Mk/elPf2q0eiEAILwRuAAAaAc7duxo9L6tgQMHNtr/rW99S2eddZY8Ho+efPJJVVRUOFEmAKCDEbgAADhN5eXlevLJJ1VTU6NRo0bpyiuvbHKMy+XSD3/4QyUlJSk/P1/PPvusA5UCADoaLz4GAAAAgCBhhgsAAAAAgoTABQAAAABBQuACAAAAgCAhcAEAAABAkBC4AAAAACBICFwAAAAAECQELgAAAAAIEgIXAAAAAAQJgQsAAAAAgoTABQAAAABBQuACAAAAgCAhcAEAAABAkBC4AAAAACBICFwAAAAAECQELgAAAAAIEgIXAAAAAAQJgQsAAAAAgoTABQAAAABBQuACAAAAgCD5/yGDdBW7obMpAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-5, 5, 201)\n",
"\n",
"plt.plot(x, np.exp(-x) - x)\n",
"plt.xlabel('x')\n",
"plt.ylabel('y')\n",
"\n",
"plt.legend([r\"$e^{-x} - x$\"])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "79126eb7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 0.5000, err_a = 1.000e+00\n"
]
}
],
"source": [
"# First Step\n",
"x0 = 0\n",
"x1 = x0 - (np.exp(-x0) - x0) / (-np.exp(-x0) - 1)\n",
"print(\"x1 = {:.4f}, err_a = {:.3e}\".format(x1, abs((x1 - x0)/x1)))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "26f21435",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 0.5663, err_a = 1.171e-01\n"
]
}
],
"source": [
"# Second step\n",
"x0 = x1\n",
"x1 = x0 - (np.exp(-x0) - x0) / (-np.exp(-x0) - 1)\n",
"print(\"x1 = {:.4f}, err_a = {:.3e}\".format(x1, abs((x1 - x0)/x1)))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c1788e25",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 0.5671, err_a = 1.467e-03\n"
]
}
],
"source": [
"# Thrid step\n",
"x0 = x1\n",
"x1 = x0 - (np.exp(-x0) - x0) / (-np.exp(-x0) - 1)\n",
"print(\"x1 = {:.4f}, err_a = {:.3e}\".format(x1, abs((x1 - x0)/x1)))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "537b6504",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 0.5671, err_a = 2.211e-07\n"
]
}
],
"source": [
"# Fourth step\n",
"x0 = x1\n",
"x1 = x0 - (np.exp(-x0) - x0) / (-np.exp(-x0) - 1)\n",
"print(\"x1 = {:.4f}, err_a = {:.3e}\".format(x1, abs((x1 - x0)/x1)))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "66589001",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" converged: True\n",
" flag: converged\n",
" function_calls: 41\n",
" iterations: 39\n",
" root: 0.5671432904109679\n",
" method: bisect"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Bisection Method\n",
"from scipy.optimize import root_scalar\n",
"\n",
"f = lambda x: np.exp(-x) - x\n",
"root_scalar(f, bracket=[0, 1], method='bisect')"
]
},
{
"cell_type": "markdown",
"id": "cfbdc1b9",
"metadata": {},
"source": [
"### Quadratic Convergence \n",
"엄밀해를 $x_r$ 로 가정하자 ($f(x_r) = 0)$. $x_i$를 중심으로 한 Taylor Expansion식을 적용하면 다음과 같다.\n",
"\n",
"$$\n",
"0 = f(x_{r}) = f(x_{i}) + f'(x_i)(x_{r} - x_i) + \\frac{f''(\\xi)}{2!} (x_{r} - x_i)^2.\n",
"$$\n",
"\n",
"이 식에서 Newton-Raphson 식을 빼면 다음과 같다.\n",
"\n",
"$$\n",
"\\begin{align}\n",
"0 &= f(x_{i}) + f'(x_i)(x_{r} - x_i) + \\frac{f''(\\xi)}{2!} (x_{r} - x_i)^2 \\\\\n",
"&- f(x_{i}) - f'(x_i)(x_{i+1} - x_i) \\\\\n",
"&= f'(x_i)(x_r - x_{i+1}) + \\frac{f''(\\xi)}{2!} (x_{r} - x_i)^2\n",
"\\end{align}\n",
"$$\n",
"\n",
"엄밀해로 부터 오차를 다음과 같이 정의하면\n",
"\n",
"$$\n",
"E_i = x_r - x_i,\n",
"$$\n",
"\n",
"위 식을 아래와 같이 정리할 수 있다.\n",
"\n",
"$$\n",
"0 = f'(x_i) E_{i+1} + \\frac{f''(\\xi)}{2!} E_i^2\n",
"$$\n",
"\n",
"해가 수렴하는 경우 $x_i, \\xi \\rightarrow x_r$ 이므로,\n",
"\n",
"$$\n",
"E_{i+1} \\approx \\frac{-f''(x_r)}{2f'(x_r)} E_i^2 = O(E_i^2).\n",
"$$\n",
"\n",
"즉 현재의 오차는 이전 단계에서의 오차에 제곱에 비례한다. \n",
"\n",
"#### 예제\n",
"위 예제에서 참값은 $x_r=0.5671432904109679$ 이다. 오차 $E_i$ 의 변화를 계산해보자."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "93ad8a06",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.5671432904109679\n",
"0.06714329041096789\n",
"0.0008322872137497273\n",
"1.253761057196101e-07\n",
"1.1868284133242923e-12\n"
]
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : np.exp(-x) - x\n",
"fp = lambda x: -np.exp(-x) - 1\n",
"\n",
"x0 = x = 0\n",
"xr = 0.5671432904109679 \n",
"\n",
"Err = []\n",
"itmax = 5\n",
"for i in range(itmax):\n",
" # Newton-Raphson\n",
" xn = x - f(x) / fp(x)\n",
" \n",
" # Compute \n",
" Ei = xr - x\n",
" Err.append(Ei)\n",
" \n",
" x = xn\n",
" print(Ei)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "51f697cd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.float64(0.1809481283173079)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-np.exp(-xr) / (2*fp(xr))"
]
},
{
"cell_type": "markdown",
"id": "b37dfde7",
"metadata": {},
"source": [
"위 오차식을 예제에 적용하면\n",
"\n",
"$$\n",
"E_{i+1} \\approx \\frac{-f''(x_r)}{2f'(x_r)} E_i^2 = 0.18095 E_i^2.\n",
"$$\n",
"\n",
"초기 오차 $E_0 = x_r - 0$ 이므로, 이 식으로 계산한 오차는 다음과 같다.\n",
"\n",
"- 매 계산마다 정확한 유효자릿수가 2배가 된다."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "5bfe9e19",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.06714329041096789 0.058202841070737574\n",
"0.0008322872137497273 0.0008157626708729342\n",
"1.253761057196101e-07 1.2534442801669386e-07\n",
"1.1868284133242923e-12 2.8443834288658173e-15\n"
]
}
],
"source": [
"# 실제 오차, 추정식\n",
"for i in range(1, 5):\n",
" print(Err[i], 0.18095*Err[i-1]**2)"
]
},
{
"cell_type": "markdown",
"id": "6b74d052",
"metadata": {},
"source": [
"### 예제 2 \n",
"아래 함수의 근을 구하시오. 초기값은 0.5로 하자.\n",
"\n",
"$$\n",
"f(x) = x^{10} - 1\n",
"$$\n",
"\n",
"Newton-Raphson 법을 적용하면 다음과 같다.\n",
"\n",
"$$\n",
"x_{i+1} = x_{i} - \\frac{f(x_i)}{f'(x_i)} = x_i - \\frac{x_i^{10} - 1}{10 x_i^9}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "79ea1fb2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 51.6500, err = 5.115e+01\n",
"x1 = 46.4850, err = -5.165e+00\n",
"x1 = 41.8365, err = -4.648e+00\n",
"x1 = 37.6529, err = -4.184e+00\n",
"x1 = 33.8876, err = -3.765e+00\n",
"x1 = 30.4988, err = -3.389e+00\n",
"x1 = 27.4489, err = -3.050e+00\n",
"x1 = 24.7040, err = -2.745e+00\n",
"x1 = 22.2336, err = -2.470e+00\n",
"x1 = 20.0103, err = -2.223e+00\n",
"x1 = 18.0092, err = -2.001e+00\n",
"x1 = 16.2083, err = -1.801e+00\n",
"x1 = 14.5875, err = -1.621e+00\n",
"x1 = 13.1287, err = -1.459e+00\n",
"x1 = 11.8159, err = -1.313e+00\n",
"x1 = 10.6343, err = -1.182e+00\n",
"x1 = 9.5708, err = -1.063e+00\n",
"x1 = 8.6138, err = -9.571e-01\n",
"x1 = 7.7524, err = -8.614e-01\n",
"x1 = 6.9771, err = -7.752e-01\n",
"x1 = 6.2794, err = -6.977e-01\n",
"x1 = 5.6515, err = -6.279e-01\n",
"x1 = 5.0863, err = -5.651e-01\n",
"x1 = 4.5777, err = -5.086e-01\n",
"x1 = 4.1199, err = -4.578e-01\n",
"x1 = 3.7079, err = -4.120e-01\n",
"x1 = 3.3371, err = -3.708e-01\n",
"x1 = 3.0034, err = -3.337e-01\n",
"x1 = 2.7031, err = -3.003e-01\n",
"x1 = 2.4328, err = -2.703e-01\n",
"x1 = 2.1896, err = -2.432e-01\n",
"x1 = 1.9707, err = -2.189e-01\n",
"x1 = 1.7738, err = -1.968e-01\n",
"x1 = 1.5970, err = -1.768e-01\n",
"x1 = 1.4388, err = -1.582e-01\n",
"x1 = 1.2987, err = -1.401e-01\n",
"x1 = 1.1784, err = -1.204e-01\n",
"x1 = 1.0833, err = -9.500e-02\n",
"x1 = 1.0237, err = -5.969e-02\n",
"x1 = 1.0023, err = -2.135e-02\n",
"x1 = 1.0000, err = -2.292e-03\n",
"x1 = 1.0000, err = -2.393e-05\n",
"x1 = 1.0000, err = -2.578e-09\n",
"x1 = 1.0000, err = 0.000e+00\n",
"x1 = 1.0000, err = 0.000e+00\n"
]
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : x**10 - 1\n",
"fp = lambda x: 10*x**9\n",
"\n",
"x0 = x = 0.5\n",
"\n",
"Err = []\n",
"itmax = 45\n",
"for i in range(itmax):\n",
" # Newton-Raphson\n",
" xn = x - f(x) / fp(x)\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "markdown",
"id": "610497dc",
"metadata": {},
"source": [
"### 예제 3\n",
"아래 함수의 해를 적절한 초기값으로 부터 근을 구하시오.\n",
"\n",
"$$\n",
"f(x) = x^2 + 10 \\sin(x)\n",
"$$\n",
"\n",
"Newton-Raphson 법을 적용하면 다음과 같다.\n",
"\n",
"$$\n",
"x_{i+1} = x_{i} - \\frac{f(x_i)}{f'(x_i)} = x_i - \\frac{x^2 + 10\\sin(x)}{2 x + 10\\cos(x)}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "348396d1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1wAAAKLCAYAAADiqYnnAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjErZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvzRIYmAAAAAlwSFlzAAAXEgAAFxIBZ5/SUgAAlddJREFUeJzs3Xd8leX9//H3dTJISALIXrINQ/YGBxsVcKFi625rraNWy9daa/ur2lpHv34dda+22moVLQ6GKCDgYq+wBBmyV1iBkH1fvz8OOSeRMALn5Drj9Xw8fHh/rjPyzpUA+eS+7+sy1lorAAAAAEDI+VwHAAAAAIBYRcMFAAAAAGFCwwUAAAAAYULDBQAAAABhQsMFAAAAAGFCwwUAAAAAYULDBQAAAABhQsMFAAAAAGFCwwUAAAAAYULDBQAAAABhQsMFAAAAAGFCwwUAAAAAYULDBQAAAABhQsMFAAAAAGFCwwUAAAAAYZLoOgCObceOHbLWOvnYdevWlSRlZ2c7+fixiDkNPeY09JjT0GNOQ485DT3mNPSY09BzNafGGDVs2PCUX0/DFcGstc4arrIZEFrMaegxp6HHnIYecxp6zGnoMaehx5yGXrTNKZcUAgAAAECY0HABAAAAQJjQcAEAAABAmNBwAQAAAECY0HABAAAAQJjQcAEAAABAmLAsPAAAQIw41eWyS18XbcttRzLmNPRCNafGmFDEOWk0XAAAAFGsuLhYeXl5KiwsPOUfRPft2ydJ8jwvlNHiGnMaeqGaU2OMkpOTlZqaqsTE8LdDNFwAAABRqri4WAcOHFBKSopq1aoln+/U7hYp/aGzuLg4lPHiGnMaeqGaU8/zlJ+frwMHDqhmzZphb7pouAAAAKJUXl6eUlJSlJaWdlrvU3qJVVVfahXLmNPQC9WcJiQkBP7M5OXlKSMj47SzHQ+LZgAAAESpwsJCpaSkuI4BRKVq1aqpsLAw7B+HhgsAACAKWWtlrT3lywiBeJeQkBD4cxRO/AkFAAAAgDCh4QIAAACAMKHhAgAAAIAwoeECAAAAgDCh4QIAAACAMKHhAgAAAIAwoeECAABA3Hr33Xd15ZVXqlOnTsrMzNSFF16o8ePHu46FGJLoOgAim7WWHdIBAEDM+uqrrzR8+HD94Q9/UM2aNTVlyhT96le/UkJCgi699FLX8RADOMOFo9g9u3Xo/Te09//9Ut5zD7uOAwAAEDbPPvusbrnlFnXt2lUtW7bUbbfdpkGDBmnixInOMmVlZem5557TzTffrB49eqhJkyZq1arVCV+Xn5+vJ554Queee65atWql7t27a+zYsdq2bdtpZ2rSpIn69Olz2u9T1pNPPqkzzzxTa9asOa332blzp1q3bq37778/RMlCizNcONre3Tr075f8x4lJsoUFMsnV3GYCAACoIgcOHFCjRo2cffynn35an376aaVek5+fr6uvvloLFixQgwYNNHz4cG3ZskXvvvuupk2bpo8//lgtWrQIT+BTsHv3br300ksaNWqUMjMzT+u9GjRooGuvvVZvvPGGfvrTn6pNmzYhShkanOHC0VpmyqRU9x8XF0nrvnWbBwAAoIqMGzdOWVlZuv76651l6NGjh37961/rn//8p5YsWXJSr3n22We1YMEC9ejRQ19++aVeeuklTZw4UX/84x+1Z88e/c///M9pZZo1a5befffd03qPsp599lnl5ubql7/8ZUje77bbbpPneXriiSdC8n6hRMOFo5jERCV37Bqo7aolzrIAAABUlU8//VT33XefHn/8cXXq1OmU3mPz5s1q0qSJfvWrX51yjjvuuEP33HOPhg0bpnr16p3w+UVFRfrHP/4hSXrkkUeUlpYWeOwXv/iF2rdvrzlz5igrK+uUM7Vp0yZkZ8jy8vL03nvvqX379jr77LND8p6NGjVS//79NWXKFO3evTsk7xkqNFyoUHLnXoFju3KpwyQAAAAn74477lCTJk30zDPPHPXYvHnz1KpVK3Xu3Fnff/99ucc++ugj3XbbbXr00Ud19dVXV1Ha0Jg3b54OHDigFi1aqGPHjkc9PnLkSEnS1KlTy41/9913uvPOO9W/f3+1atVKnTp10rBhw/THP/5RO3fuLPfciu7hKm0ur7zySuXl5emRRx5R79691bJlS51zzjl6/vnnZa09Ks+ECROUk5Oj0aNHV/j5HO9rOHfu3GN+DS+//HIVFRVp3LhxFb6vKzRcqFByl2DDpU3rZA/luAsDAABwku655x4lJibqlVde0cGDBwPja9eu1U9+8hP5fD69+eab5c7WvPXWW/r1r3+tp556KuqaLUlauXKlJFXYbEkKnK0rfZ4kLVu2TBdeeKE++OAD1a5dWxdeeKG6deumoqIivf7661q3bt1Jf/zCwkJdc801euutt3TWWWepf//+2rFjhx555BH99a9/Per506ZNkyT169evwvc73tfwxhtvrPBrWPb9pk+fftLZqwINFyqU2KyVfLVq+wtrpdXL3AYCAAA4CS1bttTVV1+t/fv367XXXpPkX6Dh+uuv18GDB/XSSy+pa9eugee//PLL+t3vfqeHHnpI/fr1065du7Rr1y7t2bPH0WdQeVu3bpWkYy70UTpe+jxJev3115Wfn69XXnlFEydO1AsvvKA333xTM2fO1MyZM9W6deuT/vgLFy6UMUZffvml3nrrLb311lsaP368EhMT9eqrryo3N7fc8+fPn6+kpCR16NChwvc71tfwxz/+sXJyco76GpZq3ry5ateurSVLlqigoOCk84cbDRcqZIwpd5aLywoBAEC0uPvuu5WSkqJXX31VO3bs0I033qhNmzbpscce09ChQ8s99+9//7tKSkp03333qVu3boH/RowY4Sh95R0+fFiSlJqaWuHj1av7F0Mr2/iUNpTnnHPOUc8/66yz1KBBg5P++D6fT//7v/+r2rVrB8a6dOmiQYMGKS8vT0uXBn+OzM7O1q5du9S0aVNVq3bsVbCP9TX861//etTXsKzWrVuroKCgUmfowo1l4XFMyZ17KX+Wf0lSFs4AACB6WGulvNwTP7H0+Yn+HwltcXG4Ip2c1DQZY077bRo3bqzrr79er776qoYNG6a9e/dq7Nixuuaaa4567ty5c0/54zz33HNau3ZtubHSpmbu3Ln61a9+Jc/zyj3+y1/+MuTLlpfeJ3WsuavoPqpOnTrp888/11133aW77rpLXbp0kc93audizjzzzArPiJXuHbZr167AWHZ2tiSpVq1ax33Pir6G99xzj6677joVH+f7tPR9I+kMJQ0Xjqla557BYvcO2d07ZOo1dBcIAACcnLxceXcd3VwcS2EYo1SG75m3perpIXmvW265Ra+99pr27t2rMWPGnPay6BWZOXOmZs+eXeFj33///VGLOkjSmDFjQt5wla5KWHqm64fy8vLKPU/yL6M+b948TZ06VVOnTlWNGjXUrVs3DR06VGPGjFF6+sl/HY51KWPpxyt7eV/pPVllsxzLD7+Gv/nNb074moyMjHIfJxLQcOGYEuo1kBo2kXb4r/e132bRcAEAgIhnrdVDDz0UOLOTmBieH3nff//9o8Y2b96svn376uqrr9bf/va3456NCZUmTZpIkrZv317h46Xjpc+T/I3Je++9p/nz52vq1KmaPXu2vvrqK82aNUvPPfecxo8ff9LLwFfmrGRpQ3To0KHjPu9Uv4Y5OTnlPk4k4B4uHJdp3yVYrOI+LgAAEPkefPBBTZw4UcOGDVPdunU1btw4rV+/3nWssCldfGL58uUVPr5smX/xs/bt25cbN8aod+/e+v3vf6+JEydq8eLFuuyyy7Rz50499thjYclat25dSdL+/fuP+7xT/RoeOHBAklSnTp3TzhoqnOHCcZn2XWVnTJYk2VVLZT1P5hSv7wUAAFUkNc1/ed5JKj17UBVnY44r9cSXmZ3Iyy+/rNdee03dunXTiy++qDfffFN/+tOf9MQTT+iFF14IQcjI06tXL9WoUUPff/+9li9fftTy8JMmTZKk4y42IfmblLFjx+rDDz/Ut99+G5asdevWVf369bVlyxbl5eVVuNDHsb6Gjz/+uF5++eXjvv/atWuVkpJSqVUWw42fnHF8bTtK5si3yaEcacv3TuMAAIATM8bIVE+Pvv9Oc8GMjz/+WH/+85/VokULvfHGG0pNTdUNN9yg+vXr6+OPPy63D1UsSU5O1k033SRJ+v3vf1/uXq6XX35Zq1atUu/evcstpf7mm29q06ZNR73XjBkzJJW//DDUevfureLi4grPyB3va/jRRx9pxYoVx3zf77//Xvv27VPXrl2PuwJiVaPhwnGZ6ulSi+CNnZbLCgEAQASaPXu27r77btWuXVtvvfVW4JKy1NRU3XHHHbLW6vHHH3ec8uRMmzZNo0aNCvwn+TcXLjtWunlwqbvuukvdunXTggULdO655+rWW2/VqFGj9Kc//UlnnHGGnnzyyXLP/9e//qV+/fpp0KBB+vnPf67bbrtNw4cP1wMPPKCUlBT9+te/DtvnN2TIEEnSN998U278ZL6Gjz766DHft3QBk8GDB4cp+amh4cIJmfZdA8csDw8AACLNmjVr9LOf/Uw+n09vvPHGUYs9XHfddWrYsKGmTZumBQsWuAlZCXv27NHixYsD/0n+RSTKjv1w2fOUlBS99957uvvuu5WamqpPP/1UW7Zs0VVXXaVPP/1ULVu2LPf83/zmN/rRj34kY4y+/vprTZ06VXl5ebr22ms1bdo09ezZU+Fy8cUXq0aNGvrwww8DYyfzNWzUqJGmTp16zK/hBx98oKSkJI0ZMyZs2U+FsRUtzI+IsH379gr3TagK9erVk+Tf1duuXibvid/7H0hOlu/p/8gkJTnJFc3KzilCgzkNPeY09JjT0GNO/ay12rNnj+rUqXPal+JFzD1cMYQ5Pb4HHnhAr732mqZMmaJOnTqd1GuON6fbtm1Tnz59NHLkSL300ksn9X4n+2fIGHPMpe9PBme4cGKt2knJyf7jwkJpfXhuogQAAEB8uPPOO5WWlqbnnnsuJO/30ksvyefz6Z577gnJ+4VSzKxSuH79emVlZWnt2rX67rvvtG/fPiUlJemtt9467utmzZqlKVOmaMuWLUpMTFRmZqZGjx6ttm3bHvM1q1ev1vjx47VmzRoVFxeradOmuuCCCzRw4MAQf1aRwSQlSWedLa04ckp75VKZtif3mwgAAADgh+rWravbbrtNTz75pNasWaPMzMxTfq+dO3fqrbfe0jXXXBPyTaVDIWYarvfff7/S1+S+8cYbmjRpkpKTk9W5c2cVFRUpKytLS5cu1dixY9W7d++jXjNv3jw9+eSTstaqffv2ysjI0PLly/XCCy9o48aNuvHGG0P1KUUU076rbGnDtWqJdPl1bgMBAAAgqv36178OyeIcDRo00Lp160KQKDxipuHKzMxUixYt1Lp1a7Vu3Vq33HLLcZ+/fPlyTZo0SRkZGXr44YcD12WuWbNGDz74oF544QV16NBB6enpgdccOnRIL7zwgjzP0//8z/+oT58+kvwbt/3xj3/UpEmT1KNHj6P2PogFpn0XBe4m+36t7OFD/hUMAQAAABxTzNzDddlll2nMmDHq0aOHatWqdcLnT5gwQZI0evTocjfBZWZmatiwYTp8+HBgH4JSn3/+uQ4fPqyePXsGmi1JqlWrlq67zn/GZ+LEiSH4bCJQ0xZSeg3/sfWk1RXvZA4AAAAgKGYarsooLCwMbLTWt2/fox4vHVu4cGG58dK6otd0795dSUlJWrZsmQoLC0Md2Tnj88m07xKoWR4eAAAAOLG4bLi2bdumoqIi1ahRI7ChWlml+xRs3Lix3HjpbtytWrU66jWJiYlq1qyZioqKtG3btjCkjgDlGi42QAYAAABOJC4bruzsbEmqsNmS/BvHpaWlKTc3V3l5eZKkw4cPKzc3V5JUu3btCl9XOl76/rGm7Bku7dgquze+9z8BAAAATiRmFs2ojPz8fElScuneUhWoVq2acnNzlZ+fr9TU1MBrSh871mvKvv+JjB079qix5ORkPfbYY5L8y2W6UrqxXOnmkkcK7W7YRCU7tkqS0rdsUPW2HVzEi0oVzilOC3Maesxp6DGnocec+llrtW/fPiUmJp72xselry+dW5w+5jT0Qj2n1lr5fD7Vq1fvtP8MHU9cnuGy1r/e3vEmtvQ5KC+5a3Cp/MKs+Q6TAAAAAJEvLlvu1NRUSVJBQcExn1O68EVKSkq5/5e+rnr16ke9pvT9yj73eJ588snjPp6dne2s8Sv9reHu3eUvG7QtghtC5y+eq127doX1NwKx5FhzilPHnIYecxp6zGnoMadB1loVFhYqISHhtN6n9IxBcXFxKGJBzGk4hHpOi4uLZa094e1Axphyq5pXVlye4Sq9VG/Pnj0VPp6fn6/c3FylpaUFmrPq1asHmqy9e/dW+LrScZeXAoZdu05SaYOVs1/atslpHAAA4llycvJJ38oAoLyCgoLj3mIUKnF5hqtx48ZKSkpSTk6O9uzZc9TiGRs2bJAkNWvWrNx48+bNtWrVKq1fv15NmzYt91hxcbE2bdqkpKQkNW7cOLyfgEMmLUNq1lrauFaSf3l406S541QAAMSn1NRUHThwQJL/XvJTPdNVekUNt1SEDnMaeqGa05KSEhUUFCg/P181a9YMRbTjisuGKzk5WR07dtTixYs1Z84cjRw5stzjc+bMkST16NGj3Hj37t21atUqzZkzR+eff365xxYtWqSioiJ169atSjpll0yHLrKlDdfKpdLQSx0nAgAgPiUmJqpmzZrKy8vTgQMHTvkHUZ/Pf9GT53mhjBfXmNPQC9WcGmOUnJysmjVrVsmiJnHZcEnSyJEjtXjxYo0fP17du3cPXJe5Zs0aTZs2TampqRo8eHC51wwZMkTjx4/XggULNHfuXPXp00eSdODAAf373/+WJI0aNapqPxEHTPuusp/811+sWS5bXCzDCjwAADiRmJiojIwMSaf+m3/uiws95jT0QjWnVb3+QMz8lLxo0SL997//LTdWXFys3//+94H6iiuuUPfu3SVJnTt31ogRIzR58mTde++96tSpk0pKSpSVlSXP83TnnXcqPT293Pulp6frtttu01NPPaUnn3xSHTp0UEZGhpYtW6bc3FxddNFF6tSpU/g/WdfatJeSkqWiQqkgX9qwRjqL5eEBAHDtVH+QLH0dC2GFDnMaetE6pzHTcOXk5Oi7774rN2atLTeWk5NT7vGbbrpJLVq00JQpU7Rs2TIlJCSoY8eOuuKKK9SuXbsKP07fvn310EMPafz48fruu+9UXFysJk2a6IILLtCgQYNC/4lFIJOU7G+6Vi2VdOQ+LhouAAAA4Cgx03ANHDhQAwcOrJLXtWvXTvfff3+lP1YsMe27ygYarqXSJdc4TgQAAABEnrhcFh6nz3ToEizWr5bNO+wuDAAAABChaLhwas5sKaX5b9CV50lrVrjNAwAAAEQgGi6cEuNL8G+CfIRdtcRdGAAAACBC0XDhlJn2XQPHpfdzAQAAAAii4cIpM+3L3Me1bZPs/j3uwgAAAAARiIYLp65eQ6lO/UBpv81yGAYAAACIPDRcOGXGmPJnuVZyWSEAAABQFg0XTk+ZhsuuWiprrcMwAAAAQGSh4cJpMe06B4v9e6SdW92FAQAAACIMDRdOi6lRS2raMlCzWiEAAAAQRMOF02baB89yWe7jAgAAAAJouHDayu7HpdXLZEtKnGUBAAAAIgkNF07fWR2khET/cV6utGmd2zwAAABAhKDhwmkzKalSq8xAbVcucRcGAAAAiCA0XAiJspcVsnAGAAAA4EfDhZAotwHyulWyBQXuwgAAAAARgoYLodHiLCkl1X9cXCytW+k2DwAAABABaLgQEiYxUcrsGKhZHh4AAACg4UIIlb2skPu4AAAAcDpscZFKXnxUdu0q11FOCw0XQqbcflyb18seynGWBQAAANHNfvBvadFsef/7O3kT3pG11nWkU0LDhdBpfKZU8wz/sbXS6mVu8wAAACAq2ZWLZT/7wF94nrR/r4wxbkOdIhouhIwxRqZd50DNfVwAAACoLHvwgLy/Px0caHSmzJifOctzumi4EFpl9+P6loYLAAAAJ89aK+8fz0gH9vkHEpPku+UemWrV3AY7DTRcCCnTPniGS7u2y2bvdBcGAAAAUcV+PklatiBQmyt/ItO0pcNEp4+GCyFlateTGjQJ1KxWCAAAgJNhN66Vff/vwYHOvWQGj3QXKERouBByZZeHFw0XAAAATsDmHpL30uNScbF/oOYZ8t30q6hdKKMsGi6EXLn9uL7NkvU8h2kAAAAQyfz3bT0tld6KYnzy/fwemYyaTnOFCg0XQq9tJ8kc+dY6eEDattFtHgAAAEQs+9kH0tJ5gdpcdq1M204OE4UWDRdCzqSlS81bB2qWhwcAAEBF7OrlsuPfDA506ilz4RXuAoUBDRfC4oeXFQIAAABl2eyd8l56zL+xsSTVriffT++W8cVWixJbnw0iRrmFM9Ysly0uchcGAAAAEcUW5Mt7/i/SoRz/QGKSfLf+Via9httgYUDDhfBo015KSvYfF+RL69e4zQMAAICIEFgkY8v3gTFz/R0yLTOdZQonGi6EhUlK9jddR7AfFwAAACTJTnpXWvhNoDbDLpWv/2CHicKLhgthY9p3DRzbb2m4AAAA4p1dPEf2o7eDAx26yVxxk7M8VYGGC2Fj2ncOFutXy+YddhcGAAAATtmtG+W9/lRwoH5j+W75jUxCgrtQVYCGC+HTrJVUPd1/7HnSmhVu8wAAAMAJeyhH3nMPSwV5/oGUVPl++Xv/dkIxjoYLYWN8CVK74Fkuu2qJuzAAAABwwhYXy3v5r1L2Tv+AMfL9/B6ZRme6DVZFaLgQVmUvK2Q/LgAAgPhirZX99wtSmZ8DzeXXy3Tu5TBV1aLhQliVXThDWzfKHtjnLAsAAACqlp30ruzX0wK16X2+zIVXOExU9Wi4EF71G0m16wVKlocHAACID94308uvSNimg8xNv5Ixxl0oB2i4EFbGmPKrFbI8PAAAQMyzq5bKvvlccKBBE/nuuN+/V2ucoeFC+JXdj2vVUllr3WUBAABAWNmtG+W9+KhUUuIfyKgp310PyKTXcBvMERouhF25M1x7s6Wd29yFAQAAQNjYfXvkPfOQVLr/anKyfHf+P5l6Dd0Gc4iGC2FnapwhNWkeqLmPCwAAIPbY/MPy/vYnaV+2f8D45Pv5b2RaZroN5hgNF6qEad8lcGy5jwsAACCm2OJieS89Lm3ZEBgzP/65TNc+DlNFBhouVImyDZe+zZL1StyFAQAAQMhYa2XfelFasTgwZoZfLt+gkQ5TRQ4aLlSNzLOlhAT/8eFcaeN6t3kAAAAQEnbSONmvpgZq0+McmStudJgostBwoUqYlOpSmet3uawQAAAg+nnffC770VvBgTbtZX72axkfbUYpZgJVptx9XCycAQAAENX8e209Gxyo31i+O34fl3ttHQ8NF6qMKbMfl75bKVtY4CwLAAAATp1/r63H2GvrJNBwoeq0zJSqpfqPi4uktavc5gEAAECl2f175P3tISkv1z+QnCzfL/8gU7+R22ARioYLVcYkJvoXzziC+7gAAACiS2Cvrb2le20Z+X5+j0yrtm6DRTAaLlSpcvdxraThAgAAiBa2pETey3+VNpfZa+tHP5fp2tdhqshHw4UqVW4/rk3rZHMPugsDAACAk2Ktlf3Py9LyRYExM/wy+QaPcpgqOtBwoWo1aS5l1PQfWyutXuY2DwAAAE7ITv1IdtaU4ED3/jJX3OQsTzSh4UKVMsawPDwAAEAUsYvnyL7/j+BAy0z52GvrpDFLqHrl7uNa4i4HAAAAjst+/528157wX5kkSXXqy/fL38skV3MbLIrQcKHKlduPa9d22d07nGUBAABAxeyeXfKee1gqLPQPpKbJ96s/ytQ4w22wKEPDhSpn6tSTGjYJ1HbVEndhAAAAcBR7OFfes3+WDuzzDyQkyHfrb2UaN3MbLArRcMGJsme5uKwQAAAgctjiYv/y71s3BsbMtbfJdOjqLlQUo+GCE+bsbsFiVZasV+IuDAAAAALsuNeklYsDtbnoCvnOG+4wUXSj4YIbbTtKCQn+48OHpI3r3OYBAACAvC8/k50xOVCbHufIXHa9w0TRj4YLTpiU6lKrtoHarlh8nGcDAAAg3Oy6b2Xffik40LyNzE/vZvn308TswRnTIXhZIQtnAAAAuGP375X34mNScbF/IKOmfLf/juXfQ4CGC86Uu/Fy3bey+YedZQEAAIhXtqhI3kuPSQf2+gcSEuS79T6Z2vXcBosRNFxwp0UbqXqa/7ikRFq93G0eAACAOGOtlf3Py9K6bwNj5kc/l8k822Gq2ELDBWeML0Fq3yVQszw8AABA1bKzpsh++VmgNucNlxlwkcNEsYeGC06VvazQrmThDAAAgKpi16yQfeeV4ECrtjI//oWMMe5CxSAaLjhVdgNk7dgqu2e3sywAAADxwu7N9t+3VXJkL9SateW77T6ZpCS3wWIQDRecMvUaSvUbBWrOcgEAAISXLS6W9+r/SgcP+AcSE/3NVq06boPFKBouOFd2eXitWuouCAAAQBywH/5bWrsqUJsf/0KmdTuHiWIbDRecK38f1xJZr8RdGAAAgBhms+bLfjo+UJt+g2TOG+4wUexLdB3AtTVr1ujjjz/W6tWrdejQIaWkpKhly5YaPny4+vbtW+FrZs2apSlTpmjLli1KTExUZmamRo8erbZt21Zx+hjRtpPk80meJ+UelDatl1qc5ToVAABATLF7dsv7+9PBgUZnylx7G4tkhFlcN1yzZ8/W008/LWutWrdurbPPPlv79u3TihUrtHz5cl166aW69tpry73mjTfe0KRJk5ScnKzOnTurqKhIWVlZWrp0qcaOHavevXs7+myil6meJrVqGzi1bVcukaHhAgAACBlbUiLvtSf8v9yWpORq8t36W5lqKW6DxYG4bbhKSkr0+uuvy1qru+++W/379w88tmbNGj300EP6+OOPNWTIEDVs2FCStHz5ck2aNEkZGRl6+OGH1ahRo8DzH3zwQb3wwgvq0KGD0tPTnXxO0cy07ypbpuHSiKvcBgIAAIghdtK48vdtXXubTONmDhPFj7i9h2vr1q3KyclRkyZNyjVbkpSZmakuXbrIWqv169cHxidMmCBJGj16dKDZKn3+sGHDdPjwYc2YMaNqPoEYY84us3DG2lWyBfnuwgAAAMQQu3aV7MR3A7XpO0i+/oMdJoovcdtwJZ3kHgOlZ6sKCwu1fPlySarw3q7SsYULF4YoYZxpcZaUmuY/LimW1ix3mwcAACAG2MO58l77P8l6/oF6DWWu+YXbUHEmbhuuBg0aqEGDBtq6dau++eabco+tWbNGS5cuVf369dWhQwdJ0rZt21RUVKQaNWqoTp2j9yho2bKlJGnjxo3hDx+DTEKC1K5ToLYrl7gLAwAAECPs2y9Je3b5C59Pvpv/Rya1uttQcSZu7+Hy+Xy6/fbb9fjjj+vpp5/WhAkT1KBBA+3bt0/ffvut2rRpozvvvFOJif4pys7OlqQKmy1JSklJUVpamnJzc5WXl6fU1NQq+1xihenQTXbxHEmSXcEGyAAAAKfDm/+l7NxZgdpc/GOZVqyqXdXituGSpPbt2+vBBx/UE088oXXr1mndunWSpNTUVHXq1ElnnHFG4Ln5+f57ipKTk4/5ftWqVVNubq7y8/NPquEaO3bsUWPJycl67LHHJEl169at1OcTSqWNZr169arsYxafO1jZb73oL7ZvVm1jlVC3fpV9/HBzMaexjjkNPeY09JjT0GNOQ485DT3Xc1qyb4+y//NyoE7q0FW1b7jNf1VRlHI9p6cqrhuur776Si+++KLOOuss3X333WratKn27dunCRMmaPz48Vq+fLkefPBBJSYmylorScfdp6D0OTg1iY2aKqFBE5Xs3CpJKlg6X9WHjHScCgAAILpYa5XzwmOyB3MkSSYlVTV/9YeobraiWdw2XNu3b9fzzz+vWrVq6b777lNKin8PgkaNGumWW27Rvn37tHDhQs2cOVNDhw4NnLEqKCg45nsWFhZKUuC9TuTJJ5887uPZ2dnOmrjS3xzs3r27Sj+u17aTdKThOjhnlnI7x86+Zq7mNJYxp6HHnIYecxp6zGnoMaeh53JOvW+my87/Kjhw5U+0LyFZivKvr6s5NcaUW6G8suJ20Yyvv/5aJSUl6tKlS4UNUr9+/SRJK1askBS8vG/Pnj0Vvl9+fr5yc3OVlpbG/Vunoezy8HbVElmvxGEaAACA6GL37pZ959XgQIduMudf4C4Q4rfh2rt3rySpevWKV2kpHT906JAkqXHjxkpKSlJOTk6FTdeGDRskSc2asYHcaWnfRSo93X3ooPT9Wrd5AAAAooS1Vt6/XpDyDvsHUqvLd+Mvj3tLDMIvbhuuWrVqSVJgoYwfWrvW/4N+6anL5ORkdezYUZI0Z86co55fOtajR49QR40rJrW61Lp9oLbLFzlMAwAAED3svC+k5cE9Yc3VP5epHV0LTMSiuG24evbsKUlatWqVPvvss3KPrVmzRpMmTZJUfpPjkSP9CziMHz9e27dvL/f8adOmKTU1VYMHs2v36TIduweO7QoaLgAAgBOxB3PKX0p4djeZ/vxcGgnidtGMVq1a6eKLL9aECRP02muv6dNPP1WTJk20b98+rVmzRtZaDR06VJ07dw68pnPnzhoxYoQmT56se++9V506dVJJSYmysrLkeZ7uvPNOpaenO/ysYoM5u7vs+Df9xYbvZA/lyKTXcBsKAAAggtn3XpcO+VclVHI1+a67nUsJI0TcNlySdP3116tt27aaOnWq1q9fr23btiklJUUdOnTQkCFDdO655x71mptuukktWrTQlClTtGzZMiUkJKhjx4664oor1K5dOwefRQw6s6VU8wzpwD7JerIrl8j0Pt91KgAAgIhkVyyWnT0jUJvLr5Op28BhIpQV1w2XJPXu3Vu9e1du6fGBAwdq4MCB4QkEGWP8Z7m+me4fWL5IouECAAA4ii0okPev54MDLc6SGTzKXSAcJW7v4UKE+8F9XNbzHIYBAACITHbye9KeXf4iIUG+G34p42OD40hCw4WIZNp3kcyRb8+c/dKW713GAQAAiDh2x1bZz8YHajP0UpkzWzpMhIrQcCEimfQaUsuzArUts8QpAABAvLPWyvvPK1JxsX/gjLoyo652GwoVouFCxDIdg3uasTw8AABAGYu+kVYuDpS+q2+WSUl1GAjHQsOFiFV2Py6t+1b2cK67MAAAABHC5ufJe+e14MDZ3aTu/dwFwnHRcCFyNW8tpWf4j0tKpG+z3OYBAACIAHbye9L+Pf4iMVG+H/+CPbciGA0XIpbxJch06BaouY8LAADEO7t7h+zUjwK1GT5apkFjh4lwIjRciGw/uI/LWuswDAAAgFv2v29IxUX+olZtmRFXug2EE6LhQkQzZ3cNFnuzpe2bnWUBAABwya5ZLrvw60BtRt8oUy3FYSKcDBouRDRT4wypWetAbZezWiEAAIg/1iuR926ZhTJaZsr0GeAuEE4aDRciXtnVCrmPCwAAxCM7e4a0aX2g9l19s4yPH+WjAV8lRLyy+3HpuxWyBfnuwgAAAFQxW5Av+8G/A7Xpfb5M63YOE6EyaLgQ+Vq1lVLT/MfFxdLqZW7zAAAAVCE77WPpwF5/kZQsM/pGt4FQKTRciHgmIUFq3yVQcx8XAACIF/Zgjuyn4wO1GTxKpk49h4lQWTRciArcxwUAAOKRnTxOyjvsL6qny1zEMvDRhoYLUcGcHWy4tHuH7K5t7sIAAABUAbt7h+yMyYHajLhKJi3dYSKcChouRAVTu67UpHmgtss4ywUAAGKb/egtqaTYX9SuKzN4pNtAOCU0XIgaZVcrtFkLHCYBAAAIL7t5g+zcWYHaXHqtTFKyw0Q4VTRciBqmc69gsWaZbP5hd2EAAADCyPvorWDRpLlM34HOsuD00HAherRuJ1U/ct1ycbG0cqnbPAAAAGFgN3wnLZ0XqH2XXivjS3CYCKeDhgtRwyQk/OCywnnHeTYAAEB08j4uc3areRupax93YXDaaLgQXTr3DBzarAWynucwDAAAQGjZtaukMnuO+i69RsYYh4lwumi4EFVMxx6S78i37cED0sa1bgMBAACEULl7t1q1lcpc3YPoRMOFqGLS0qU2HQK1zZrvMA0AAEDo2NXLpG+zArXv0ms5uxUDaLgQdcquVkjDBQAAYoU34Z1gkXm21L6LuzAIGRouRJ1yy8NvWi+7N9tdGAAAgBCw362UVi8L1L5LOLsVK2i4EH0aNpHqNQyUdhmbIAMAgOjmTXo3WJzVQaZtR3dhEFI0XIg6xhguKwQAADHDbvhOWrE4UPtGXe0wDUKNhgtRyXTpHSxWLZUtKHAXBgAA4DR4k8cFi5aZUvuuzrIg9Gi4EJ3O6iClpPqPiwql1VnHfz4AAEAEsls2SEvmBmrfyDHcuxVjaLgQlUxiknR2t0Btl3JZIQAAiD520nvBomlLqeziYIgJNFyIWj+8j8ta6zANAABA5dgdW2UXfh2ofaM4uxWLaLgQtUzHHlLpX0r790ibN7gNBAAAUAl26odS6S+MGzaRuvVzmgfhQcOFqGVq1JJatQ3UrFYIAACihT2wT/abzwO1uWC0jI8fzWMRX1VENdOpZ+CYhgsAAEQLO32CVFzkL2rWlukz0GkehA8NF6Ka6VLmxtLvv5PN2ecuDAAAwEmw+YdlZ34SqM3Qi2WSkhwmQjjRcCG6NWkh1a7rP7ZWdtlCp3EAAABOxH7xmZSX6y9Sq8ucf6HbQAgrGi5ENWPMUasVAgAARCpbXCQ79aNAbc6/UKZ6msNECDcaLkQ907l3sFixRLaoyF0YAACA47DzvvSvrixJiYkyQy92GwhhR8OF6Neuk5RczX9ckCetznKbBwAAoALW2vJnt/oMlKlVx2EiVAUaLkQ9k5QsdegWqO2SuQ7TAAAAHMPqZdKW4L6hZthl7rKgytBwISaYbn0Cx3bJPFnPc5gGAADgaN60j4NFh24yTZq5C4MqQ8OFmGA695LMkW/nA3ul779zGwgAAKAMu3ObVGZxL9+wSxymQVWi4UJMMOk1pMyzA7VdMsdhGgAAgPLs9I8la/1FozOls7u7DYQqQ8OFmGG6lrmscDH3cQEAgMhgcw/Jfj09UJuhF8sY4zARqhINF2JG2YZLO7bIbt/iLgwAAMAR9stPpcICf5GeIdN3kNtAqFI0XIgZpm4D6cyWgZrLCgEAgGu2pER2xqRAbc6/SKZ0OxvEBRouxBTTtW/g2C6m4QIAAI4tnSvtzfYfJyTIDLrIbR5UORouxBTTLdhwacMa2dKd3AEAABzwPi9zdqt7fzY6jkM0XIgtTVtIdeoHSjZBBgAArhRtWu/f7PgIM3iUwzRwhYYLMcUYU+4sF5cVAgAAVw5P/m+waNZKat3OXRg4Q8OFmFPussJvs2QP5bgLAwAA4pKXe0j5Mz8J1GbQSJaCj1M0XIg9bdpLNWr5jz1Pduk8p3EAAED8yft8kmx+nr9Iy5Dpfb7bQHCGhgsxx/gSyl9WuPAbh2kAAEC8sZ6nw58ELyc05w5jKfg4RsOFmGS69w8Wq5bI5h12FwYAAMSXVUtVsm2z/9gYmYEsBR/PaLgQmzI7SukZ/uPiYtms+W7zAACAuOHNCt67pc69ZOo2cBcGztFwISaZxESZLn0CtV3EZYUAACD87L49Upn7x30DOLsV72i4ELNMjzKXFS5fKFuQ7y4MAACIC/arqZLnSZIS6jeSzu7qNhCco+FC7GrXRUpN8x8XFkrLF7rNAwAAYpotKZH98rNAnTr8UhlfgsNEiAQ0XIhZJilJpkuvQM1qhQAAIKyWLZD2ZfuPExKUOmSU2zyICDRciGllVyu0WQtkCwscpgEAALHMmzUlcJzSZ4ASzqjjMA0iBQ0XYtvZ3aRqqf7jgjxp+SK3eQAAQEyyu3dIK4I/Z6ReeLnDNIgkNFyIaSa5mkyX3oHaLvjKYRoAABCr7JefSdb6iwZNlNyph9tAiBg0XIh5pte5gWO7dB6rFQIAgJCyJSWy33weqM35w2WMcZgIkYSGC7Hv7O5lVisskM1a4DYPAACILcsXSQf2+o8TEmX6DXabBxGFhgsxzyQlyXQtswnygi8dpgEAALHG+2pqsOjSWyajprswiDg0XIgLptd5wWLZQtn8w+7CAACAmGFz9knL5gdq37nDHKZBJKLhQnxo30VKy/AfFxXKLpnnNg8AAIgJdvZMqaTEX9SqI53d1WUcRCAaLsQFk5go071foGa1QgAAcLqstbJfTwvUpv8QGV+Cw0SIRDRciBumZ3C1Qi1fJHv4kLswAAAg+q1fLW3fHCjNOUMchkGkouFC/GjbSSq9ibWkWHbxHLd5AABAVLNlF8to20mmfiN3YRCxaLgQN0xCgkyP/oHazp3lMA0AAIhmNj9Pdn7wFgVz7lCHaRDJEl0HiAT79+/Xhx9+qMWLFys7O1vJycmqX7++OnXqpOuuu+6o58+aNUtTpkzRli1blJiYqMzMTI0ePVpt27Z1kB6VYfoMkJ35ib/4Nkt2/x6ZWnXchgIAAFHHLvxaKsjzF6lpMt37H/8FiFtxf4ZrzZo1+vWvf63JkycrISFBPXv2VGZmpg4dOqSJEyce9fw33nhDzz//vDZv3qxOnTqpTZs2ysrK0gMPPKB581j5LuK1bi/Vqe8/tlZ2HntyAQCAyrNflVkso/d5MsnVHKZBJIvrM1x79+7Vo48+qqKiIt1zzz3q3bt3ucfXrl1brl6+fLkmTZqkjIwMPfzww2rUyH+d7po1a/Tggw/qhRdeUIcOHZSenl5lnwMqxxjjP8s1+T1JRy4rHH6Z21AAACCq2B1bpLUrA7Vh7y0cR1yf4Xr77beVm5ur66677qhmS5LatGlTrp4wYYIkafTo0YFmS5IyMzM1bNgwHT58WDNmzAhvaJw202dAsNi0Tnb7FndhAABA1Cl7dktNW0jN2xzzuUDcNlyHDh3S7NmzVb16dQ0ePPiEzy8sLNTy5cslSX379j3q8dKxhQsXhjYoQs40biY1axWo7dyZ7sIAAICoYouLZWd/HqjNOUNljHGYCJEubi8pXL16tYqKitSpUyclJiZqzpw5+vbbb1VcXKwmTZqoX79+qlWrVuD527ZtU1FRkWrUqKE6dY5eZKFly5aSpI0bN1bVp4DTYPoMkN20XpL/skJ76bX8ZQkAAE5s+UIpZ7//ODFRpu9Al2kQBeK24dq82b9JXc2aNfXHP/5Ra9asKff422+/rdtvv139+vWTJGVnZ0tShc2WJKWkpCgtLU25ubnKy8tTampqGNPjdJle58u+/0/JWil7p7TuW6lNe9exAABAhPO+LrNYRte+Muk1HKZBNIjbhis3N1eS9MUXXygxMVG33nqrevbsqfz8fE2ZMkUTJ07Us88+q8aNG6t58+bKz8+XJCUnJx/zPatVq6bc3Fzl5+efVMM1duzYo8aSk5P12GOPSZLq1q17Kp9aSCQm+r816tWr5yxDWNWrp70du6twmf8S0JSsuarR7/ywfsiYn1MHmNPQY05DjzkNPeY09JjTk1Oyb492Zy0I1LVGXqFqx5gz5jT0onVO4/YeLs/zJEklJSW68cYbNXjwYNWoUUP169fXDTfcoL59+6q4uFgfffSRJMlaK0nHveys9DmIDikDLwwc5301Tbao0GEaAAAQ6fJnTZG8EkmSr24DJXfu5TgRokHcnuEqPQNljNGAAQOOenzQoEGaM2eOVq5cWe75BQUFx3zPwkL/D+wpKSknleHJJ5887uPZ2dnOmrjS3xzs3r3bycevCvasTlJyslRYKHswR7unfyLTI3ybFsbDnFY15jT0mNPQY05DjzkNPeb05JRMnRA4tn0GKHvv3mM+lzkNPVdzaowpt0J5ZcXtGa7SL1itWrWUlJR0zMcPHDggKXh53549eyp8v/z8fOXm5iotLY37t6KESa1ebld475vpDtMAAIBIZjdvkLYGF0cz/U68yjUgxXHDVbqqYG5uboVnkQ4dOiQpeLaqcePGSkpKUk5OToVN14YNGyRJzZo1C1dkhIHpPyRYLF8oe2CfuzAAACBi2Tll9lptmSnToLG7MIgqcdtwNWvWTPXr11dhYaG+++67ox5fsWKFJKlVK/9+TcnJyerYsaMkac6cOUc9v3SsR48e4YqMcGjbSap9ZHESz5OdM9NpHAAAEHmsVyI794tAbfoNcpgG0SZuGy5JuvTSSyVJ//jHP5STkxMYX79+vSZOnChJGjZsWGB85MiRkqTx48dr+/btgfE1a9Zo2rRpSk1NPalNlBE5jM9X7pIA+810Fj8BAADlrcqSDhy5XyshQabneW7zIKrE7aIZkjRkyBAtW7ZMc+bM0d13363MzEwVFBRo9erVKi4u1pAhQ9S3b9/A8zt37qwRI0Zo8uTJuvfee9WpUyeVlJQoKytLnufpzjvvVHp6usPPCKfC9B8iO2mcv9i2Sdq4VmpxlttQAAAgYpS7nLBjD5kM9t7CyYvrhsvn8+nuu+/W1KlT9fnnnwcuI2zdurWGDRum888/el+mm266SS1atNCUKVO0bNkyJSQkqGPHjrriiivUrl27qv4UEAKmfiPprA7Sd/4VKe3X02VouAAAgCSbnye7aHag9nE5ISoprhsuyd90XXDBBbrgggtO+jUDBw7UwIEDwxcKVc70HyJb2nDNmyU75qcyScfe5BoAAMQHu3iOVHhkW6DUNIm9t1BJcX0PF1DK9DxHSq7mLw7n+v9yBQAAca/s5YSm5zn8QhaVRsMFSDIp1WV6nhuo7RefOkwDAAAigd2/x79gxhGmL5cTovJouIAjzPllLitdvUx2x1Z3YQAAgHN27heS9fxFnfpSm/ZuAyEq0XABpVq1lZo0D5T2q88chgEAAK6Vu5yw70AZHz86o/L4rgGOMMbInBc8y2W/ni5bVOQwEQAAcMVu2SBt+T5QczkhThUNF1CG6TtQKr0Z9lCO7JK5TvMAAAA37OyZwaJlpkzDJs6yILrRcAFlmLR0mR7nBGr7JYtnAAAQb6xXIjtvVqA27L2F00DDBfxAucUzVi2V3bXdXRgAAFD1vs2S9u/1HyckyPQ8z20eRDUaLuCH2rSXGp0ZKFkiHgCA+GJnBxfLUMceMhk13IVB1KPhAn7AGCNz/vBAbb+eKlu6wzwAAIhpNj9PdtHsQO3jckKcJhouoAKm/xApuZq/OHRQdv5XbgMBAIAqYZfMlUp/0ZqaJnXu5TYQoh4NF1ABUz293PKv9vOJstY6TAQAAKqCnfdF4Nj06C9TunoxcIpouIBjMINGBItN66T1q92FAQAAYWcP5UgrFwdq0/t8h2kQK2i4gGMwTVtIbTsFavv5RHdhAABA2NlF30glJf6iZm2pbUe3gRATaLiA4/ANGhk4tgu/li1dIhYAAMQcO7fM5YQ9z5HxJThMg1hBwwUcT9c+0hl1/cclJbJffuY2DwAACAu7b4/03YpAzeWECBUaLuA4TEKCzIALA7Wd9YlscZHDRAAAIBzs/C+l0gWy6jWUWma6DYSYQcMFnIA5/wIpMclfHNhXbvUiAAAQG8qtTtjrPBljHKZBLKHhAk7AZNSUKbPpof3sQ5aIBwAghtid26SNawO16T3AYRrEGhou4CSYYZcGi60bpRWL3IUBAAAhVe7qlSbNZZo0cxcGMYeGCzgJptGZ5Xaa9z770F0YAAAQMtba8pcTslgGQoyGCzhJvuGXB4tVS2U3rXMXBgAAhMbmDdKOLYHS9DrPYRjEIhou4GRlni21OCtQWs5yAQAQ9cpdTtiqrUy9hu7CICbRcAEnyRgjU+Ysl53/peye3Q4TAQCA02E9z78c/BEsloFwoOECKsF07yfVqe8vPE/2sw/cBgIAAKdu3bfS3iO/PDU+mV7nuM2DmETDBVSCSUiQGX5ZoLZffiZ7YJ+7QAAA4JSVu5ywXSeZGme4C4OYRcMFVJI5d5hU88hfyEWFnOUCACAK2ZIS2YVfB2pWJ0S40HABlWSSq5W/l2vmJ7IHDzhMBAAAKm3VUqn03+/ERP9tA0AY0HABp8AMuFBKr+EvCgtkp37kNhAAAKiUcpcTduwhUz3dXRjENBou4BSYainl7+WaMUk296C7QAAA4KTZokLZxbMDNasTIpxouIBTZAaNkEp/G5afJzttgttAAADg5CxfJOXn+Y+rpch07uU2D2IaDRdwikxKdZmhlwRqO/1j2YM5DhMBAICTUW7vrS69ZapVc5gGsY6GCzgNZsio4FmuvMOyU953GwgAAByXLSiQzZofqE3Pcx2mQTyg4QJOg6meLnPRFYHafj5Jdm+2w0QAAOC4li+QCvL9xympUsfubvMg5tFwAafJDBol1aztL4qLZCe+4zYQAAA4Jjv/q8Cx6dpHJinZYRrEAxou4DSZatVkRl0dqO3X02R3bHWYCAAAVMQW5Msu43JCVK2QN1x33nmnPvzwQx04wEawiB/m3GFSvYb+wvNkP3rLbSAAAHAUmzVfKiz0F6lpUodubgMhLoS84dq1a5f+85//6LbbbtOTTz6prKysUH8IIOKYxESZS68N1HbBV7Ib1jhMBAAAfsgu+OHlhEkO0yBehLzhuvzyy1W7dm2VlJRo7ty5+stf/qI777xTH330EWe9ENNMr/OkM1sGau+dV2WtdZgIAACUsvmHpWULA7XpdZ7DNIgnIW+4fvSjH+n555/Xvffeq27dusnn82nXrl16++23ddttt+mpp57irBdikvH55Bvzs+DA+tWyc2e5CwQAAALs0vlS0ZHLCaunS+07uw2EuJEYjjf1+Xzq0aOHevToob179+rzzz/XjBkzlJ2drTlz5mjOnDmqX7++hg4dqoEDB6pmzZrhiAFUOdOus9S9n7RotiTJ/vcN2W59ZaqlOE4GAEB8K7fZcfd+MolcToiqEfZVCmvXrq0rr7xSzz33nH73u9+pV69eSkhI4KwXYpbvyp9IiUd+l7F/j+yU/7oNBABAnLOHc6UViwI1qxOiKoXlDFdFjDHq2rWrunbtqr179+pvf/ubVq1apZKSksBZr4YNG+riiy/W4MGD5fOxYj2ik6nXUGbYZbKfvC9Jsp9+IHvuMKlePcfJAACIT3bJXKm42F+kZ0jtuJwQVadKu5rs7GyNGzdOv//977Vq1arAeIsWLeTz+bRjxw69+uqr+v3vf6+cnJyqjAaElBlxpVTzDH9RVCjvnVfdBgIAII6VW52we3+ZhASHaRBvwn6Gy/M8LViwQNOnT1dWVpY8z5Mkpaena+DAgRo2bJgaNmyo/fv367PPPtOkSZO0fv16vf3227r11lvDHQ8IC5NSXWb0jbL/eNo/sGSu8ufMVErfgS5jAQAQd2zuIWnlkkDN5YSoamFruHbt2qXp06dr5syZ2r9/f2A8MzNTw4cPV79+/ZSYGPzwtWrV0pgxY9SjRw/df//9Wrx4cbiiAVXC9Bsk+810afUySVLOK/+n5E49HacCACC+2CVzpJIjlxNm1JQyO7oNhLgT8oZrzpw5mjZtmpYvXx7Ygyg1NVXnnnuuhg8frmbNmh339a1bt1atWrXKNWlANDLGyHfd7fIe+pVUXCRvb7YOvfWSdPmNrqMBABA3yq1O2OMcLidElQt5w/XUU08Fjlu0aKFhw4bp3HPPVUrKyS+LXfbMFxDNTMMmMiOvkv3obUnS4U/Gy9e5j0zrdo6TAQAQ++yhHGnV0kDN5YRwIeSdTVJSkvr376/hw4erTZs2p/Qezz//fIhTAe6YC6+QnfeltH2zZK28N5+T7w9PyiQlu44GAEBMs4tmS0fWD1DNM6Sz2rsNhLgU8lUKX375Zd1+++2n3GwBscYkJsl3/R3BgW2bZD/8t7tAAADEiXKrE/Y4R8bH5YSoeiFvuNLS0kL9lkDUM2d1UPVRYwK1nfqR7Lds9g0AQLjYgwekb5cFai4nhCvsLgxUkYzrb1NC0xb+wlp5/3ha9vAhp5kAAIhVduE3kj1yOWGtOhL3T8MRGi6giphqKar16wel0tWR9mbL/ucVp5kAAIhV5S4n7HmOjI8fe+EG33lAFUpq3Vbm4h8HajtnprzZMxwmAgAg9tgD+6Q1KwI1lxPCJRouoIqZC68od1mD/ffzslu+dxcIAIAYYxeVuZywdj2pVVu3gRDXaLiAKmYSEuT7+T1SWoZ/oLBQ3ouPyR7OdRsMAIAYUW6z457nyhjjMA3iHQ0X4ICpU1++m8dKpf8A7Nom75/PyFrrNhgAAFHO7tsjrV0VqLmcEK7RcAGOmI49ZEb9KDiweI7sJ++7CwQAQAywC7+WSn+BWbeB1IK9YeEWDRfgkBl1tdSxe6C2H/yr3KpKAACgcsqvTsjlhHCPhgtwyPh88v1srFSvYWDMe/0p2TKXQgAAgJNj9+6W1n0bqE0vLieEezRcgGMmvYZ8v/qjVD3dP1BcJO/5h2V3bXMbDACAKGMXfB0s6jeSzmzlLgxwBA0XEAFMw6by3XG/lJjoHzh0UN4zD8nu3+s2GAAAUYTLCRGJaLiACGEyO8rcdFdwYNd2eU/+P9mc/c4yAQAQLWz2TmnDmkDN5YSIFDRcQATx9RkgM/rG4MD2zf6m61COu1AAAESBcotONWwiNWnhLAtQFg0XEGF8F10hc3GZ5eK3bpT31B9lcw+6CwUAQIQre/+W6XkelxMiYtBwARHIXPxjmYuuDA5sWi/v8fv8qy8BAIBy7K5t0sa1gZrNjhFJaLiACGSMkbn8epnhlwcHt2+W99hvZbduchcMAIAIVG51wsbNZJo0cxcG+AEaLiBCGWNkrrxJZlSZywv3Zcv7632y3610FwwAgAjzw9UJgUhCwwVEMGOMfJdeI3PtbVLpteiHD8n7vz/Im/mJrLVuAwIA4JjdsVXavCFQ03Ah0tBwAVHAN/Ai+W79rZSY5B8oKZZ960XZN/4mW1ToNhwAAA6VW52waQuZRk3dhQEqQMMFRAnTvb98//OwVLN2YMx+Pd2/mMaOLQ6TAQDgjp3/ZeCYs1uIRDRcQBQxbdrL94cnpTbtg4Mb18r7893yPp8o63nuwgEAUMXs1k3StuBiUmx2jEiU6DpAJDl06JDuvvtu5eTkqHHjxnr66aeP+dxZs2ZpypQp2rJlixITE5WZmanRo0erbdu2VRcYccnUqi3f/zwsO+512RmT/YOFhbL/eUV2yVz5brxTpk59tyEBAKgCdkHw7Jaat5Gp39hdGOAYOMNVxhtvvKGDB0+8uewbb7yh559/Xps3b1anTp3Upk0bZWVl6YEHHtC8efOqICninUlMku+aW+W7/X4po2bwgVVL5f3xdnkT3+XeLgBATLPWys4vszohZ7cQoTjDdcSyZcs0a9YsDR06VNOmTTvm85YvX65JkyYpIyNDDz/8sBo1aiRJWrNmjR588EG98MIL6tChg9LT06sqOuKY6dZXvtbt5P3reWnJXP9gYaHsR2/Jzv5cvjE/kzr3kild4RAAgFixeYO0c2ug5P4tRCrOcEkqLCzUq6++qqZNm+riiy8+7nMnTJggSRo9enSg2ZKkzMxMDRs2TIcPH9aMGTPCmhcoy9SoJd/t98v87Nflz3bt2i7vuYflPf5b2VVL3QUEACAMyl1O2Kotl9MjYtFwSXrvvfe0c+dO3XzzzUpISDjm8woLC7V8+XJJUt++fY96vHRs4cKF4QkKHIMxRr6+g+R7+EWZIRdLpswf7XXfynvy/6nk//4gu3Ixe3cBAKIelxMimsT9JYUbN27UxIkTNXDgQHXo0EG7du065nO3bdumoqIi1ahRQ3Xq1Dnq8ZYtWwbeE3DBVE+X+dHPZc8dKu+/b0rLyzT/32bJ+zbLv0fJsMtkep8nU7qvFwAA0eT776Tsnf5jY2R60HAhcsX1GS7P8/Tyyy+revXquu666074/OzsbEmqsNmSpJSUFKWlpSk3N1d5eXkhzQpUhmnaUgl3PSDfvY9JmR3LP7jle9l/PC3v3p/Ke/8fsju2VvwmAABEqLJ7b+msDjJnVPyzGRAJ4voM15QpU7R27VrdfvvtysjIOOHz8/PzJUnJycnHfE61atWUm5ur/Px8paamHvf9xo4de9RYcnKyHnvsMUlS3bp1T5gpXBIT/d8a9erVc5Yh1jiZ03oDZPudr8Lli3T4w7dVsPCb4GMHD8h++oHspx8o6eyuSj13mFL6D5Kv5hlVl+808X0aesxp6DGnocechl40zan1PO1eNFulF8jXGHiRqkdg7mia02gRrXMatw1Xdna23nnnHXXo0EEDBw48qdeU3vtyvBXfuD8GkcYYo2qdeqhapx4q2rRehz96W3lfTJXKLBtftGKJilYsUc6rTyq5Uw+lnDdUKX0HyJdew2FyAACOVrR6mbw9R24B8flUrf8gt4GAE4jbhuu1115TcXGxbr755pN+TekZq4KCgmM+p7DQ/0NsSkrKCd/vySefPO7j2dnZzhq40t8c7N6928nHj0URMaepGdKPfiHfxdfIzpkp+8Wn0rZNwce9EhUunafCpfOU8+JfpbO7yfQ4R6ZLL5m0E58FrmoRMacxhjkNPeY09JjT0IumOfWmTgwWbTtpb1GJFIG5o2lOo4WrOTXGlFudvLLituFatGiR0tLS9Nprr5UbLyoqkuRvdh588EFJ0n333aeUlJTAJX579uyp8D3z8/OVm5urtLS0E15OCLhk0jJkhlwsO3iU9P13svO+9C+vu39v8EklxVLWfNms+bI+n9S2k0zXPjJd+8rUdne5KwAgflmvRHbh14GavbcQDeK24ZKk3NxcrVy5ssLHCgsLA4+VlJRIkho3bqykpCTl5ORoz549Ry2esWHDBklSs2bNwpgaCB1jjNQyU6ZlpuxVP5HWrpKd/6X/H7ODB4JP9Dxp1VLZVUtl//OK1LyNTLe+Mt37SQ2bsrEyAKBqrFkhHdjnP05I8P87BES4uG24xo0bV+H4rl279Mtf/lKNGzfW008/Xe6x5ORkdezYUYsXL9acOXM0cuTIco/PmTNHktSjR4+wZAbCyfh8UubZMplny/7o59LqZbILvpJdMrd88yVJG9fKblwr++G/pQZNjjRf/aUWbWi+AABhU3bvLbXvKsO9xogCcdtwnaqRI0dq8eLFGj9+vLp37x64nnPNmjWaNm2aUlNTNXjwYMcpgdNjEhKkDl1lOnSVve42af1q2cVzZRfPlnbvKP/knVtlp/xXdsp//c1X34EyfQbI1GvoJjwAICbZkhLZRcHVdtnsGNGChquSOnfurBEjRmjy5Mm699571alTJ5WUlCgrK0ue5+nOO+9Uenq665hAyBhfgtSmg0ybDrJX3iRt2yS7eLbs4rnSpnXln7xzq+xHb8l+9Jb/NX0HyvQ8VyaNPxMAgNP0bZZ0KMd/nJgo07Wv2zzASaLhOgU33XSTWrRooSlTpmjZsmVKSEhQx44ddcUVV6hdu3au4wFhY4yRmjSXadJcGvUj2T27ZJfM9d/z9d0P7odcu1J27UrZd17xL7QxaKR/c0ouOQQAnIJymx2f3V2mepq7MEAl0HD9QP369Y95f1dZAwcOPOn9u4BYZerUlxlysTTkYtnsnbJzZ8nOmSnt2BJ8UnGx/16wBV9JTVvKDB4p03eQTFKSs9wAgOhii4v8l7UfYXqd5zANUDk0XABCwtRtIDNyjOyIq6RN6/z7fM2dVX7BjS0bZN98Tvbjt2WGXy5z/gUy1U68Zx0AIM6tWCIdzvUfJyfLdOntNA5QGTRcAELKGONfNr55G9krfyJlzZc3Y5K0amnwSfv3yo57XXbyezIXjpYZPEomKdldaABARLMLylxO2KmnTAr7nSJ60HABCBuTkCB166uEbn1lt22SnT5R9ptpUnGx/wmHcmTf/6fs55NkLr9OpvcA//L0AAAcYYsK/VuUHOHjckJEGX6yAVAlTONm8l1/u3yPvioz/DKp7KWEe3fLvv6UvL+MlV2/2llGAEAEWrZQys/zH1dLkTr2dJsHqCQaLgBVytSqI99VP5Xv0ddkhl4iJZQ50b5pvbzH7pX31kuypdfqAwDiml0Q3OzYdOkjU62awzRA5dFwAXDCZNSQ7+qb5fvzC+VXm7JWduZkeX+8XXbJHHcBAQDO2YJ82aXzAjWbHSMa0XABcMrUayjfLb+R7zePSI3ODD5wYJ+85x+R9+ZzsqWXkgAA4orNWiAVFviL1DTp7O5uAwGngIYLQEQwmR3l++PTMpddJ5VZsdB++Zm8P98tu2GNw3QAABfs/C8Cx6ZrH/ZwRFSi4QIQMUxiknwjx8j3x6el5m2CD+zaLu/x++TN/ETWWmf5AABVx+Yd9i+YcQSbHSNa0XABiDimYVP57vurzIgxkjny11RJsexbL8q+8axsUaHbgACAsLNL50rFRf4iLUNq38VtIOAU0XABiEgmMVG+y6+T756/SDXPCIzbr6fJe/w+2X17HKYDAISbnV9mdcLu/WQS2T4W0YmGC0BEM5lny/eHp6TW7YKDG9fKe+w3Kt68wV0wAEDY2EM50opFgZrLCRHNaLgARDxTq7Z89/xFZuCI4ODebO257xcqXLHEVSwAQJjYBV9LJSX+ouYZUtuObgMBp4GGC0BUMIlJ8l17q8yPb5GMkSTZ3IPa++Bdsgu/cZwOABBKdt6swLHpdZ6ML8FhGuD00HABiCq+waPku/W3UuKRpYGLCuW98ld58744/gsBAFHB7tklfbcyUJs+AxymAU4fDReAqGO695dv7J9l0jP8A54n+9qT8ubMdJoLAHD6bNlfoDVoUn6bECAK0XABiErmrA6q/fALMjVq+QesJ/v3p+R987nTXACA02PnlrmcsPf5MkcuIweiFQ0XgKiV1KKNav/5OSmjpn/AWtl/PiNv/pdugwEATond8r20dWOgNn25nBDRj4YLQFRLat7av1dX4EyXlX39KdkVi53mAgBUXtmzW2qZKVO/sbswQIjQcAGIeqZxM3/TVXpPV0mxvBcflV2/2m0wAMBJs55X7v4t0/t8h2mA0KHhAhATTKMz5fvVA1K1FP9AQb68v/1Jdvtmt8EAACdn7Spp727/sfGx2TFiBg0XgJhhWmbKd/v9UkKifyD3oL/pOnjAbTAAwAmVu5ywfReZmme4CwOEEA0XgJhiOnSV7+axgc2Rlb1T3guPyhYVuQ0GADgmW1wku/DrQG36cDkhYgcNF4CYY3qeK3PFjcGBtStl//W8rLXuQgEAjm3FYin3oP84KVmmWz+3eYAQouECEJPM8MtlzhkSqO3sz2WnjHeYCABwLOX23urSWya1usM0QGjRcAGIScYYmetulzLPDozZD96UXb7IYSoAwA/Z/MOyS+cGai4nRKyh4QIQs0xikny3/k6q19A/YK281/9Pds9ut8EAAAF28VypsNBfVE+XOvZwGwgIMRouADHNZNSQ77bfSUnJ/oFDB+W9/LhsMYtoAEAksHNnBo5Nz3NkEpPchQHCgIYLQMwzZ7aUue624MCGNbLj/u4uEABAkmRz9kmrlgZq02eAwzRAeNBwAYgLvv5DZM4bHqjtjEny5n/pMBEAwM7/WvI8f1G7rtSmg9tAQBjQcAGIG+bHt0jNWgdq+68XZPfscpgIAOKbnTMjcGx6nS/j40dTxB6+qwHEDZOULN+tv5VSUv0DebnyXn9S1itxGwwA4pDdukn6/rtAbfoNcpgGCB8aLgBxxdRrKHPtrcGB71bKTn7fXSAAiFN29vRg0byNTJPm7sIAYUTDBSDu+PoOkukdvDHbTviP7PrVDhMBQHyxJSWyc2YGatN/sLswQJjRcAGIS+baW6U69f2F58l77f9kCwrchgKAeLFysXRgn/84IVGmN5sdI3bRcAGIS6Z6mnw3j5XMkb8Gd++Q/eBNt6EAIE7Ybz4PFl16yaTXcBcGCDMaLgBxy7TpIHPB5YHafj5Rds0Kh4kAIPbZ3EOyS+YEal//IQ7TAOFHwwUgrplLfiw1OtNfWCvvjb9xaSEAhJGd/4VUXOwvMmpKZ3d3GwgIMxouAHHNJCXLd9OvgpcW7tou++G/3IYCgBhW9nJC02egTGKiwzRA+NFwAYh7plVbmeGXBWo7fYLs2lXuAgFAjLLbN0sb1gRqcw6rEyL20XABgCRz6TVSw6b+wlp5/35BtvSSFwBASJRbLKNZK5mmLd2FAaoIDRcA6MilhTf+MjiwdaPs1I/cBQKAGGO9Etk5MwK1YbEMxAkaLgA4wrTpIHP+BYHaTvyP7O4dDhMBQAxZuUTav9d/nJDA3luIGzRcAFCGGX2jf9UsSSoslPfWi7LWug0FADGg3OWEnXrJlP5dC8Q4Gi4AKMOkpctcfXNwYMVi2QVfuQsEADHAHj4ku7jM3lssloE4QsMFAD9gep8vdegWqO2412Xz8xwmAoDoZud/JRUX+Yv0GlLHHm4DAVWIhgsAfsAYI9+1t0qJSf6B/XtlJ41zGwoAopj9Znrg2PQZIFP69ysQB2i4AKACpn4jmQsuD9R26keyO7Y6TAQA0clu3yKtXx2oWZ0Q8YaGCwCOwVx0lVS7rr8oKZb37mssoAEAlWS/nhYsmraUadbKXRjAARouADgGU62afFf9NDiwfKGUNd9dIACIMra4qPzlhOdwdgvxh4YLAI6nxzlSu86B0nv3NdmiIoeBACCKLJ0vHTzgP05Mkuk3yG0ewAEaLgA4DmOMfD+6RfId+ety9w7Zzye6DQUAUcL74tPAseneXyYtw2EawA0aLgA4AdOkmcygkYHaThonezDHYSIAiHw2e6e0akmgNucPdxcGcIiGCwBOghl1tVQ9zV/k5cpO+I/bQAAQ4exXU6XShYbqN5YyO7oNBDhCwwUAJ8Gk15AZeXWgtrM+8S91DAA4ii0pKbc6oTl/uIwxDhMB7tBwAcBJMoNGSvUa+gvPk/f+P9wGAoBItXyhtH+v/zghUabfYLd5AIdouADgJJmkJPmuuCk4kDVfdtVSZ3kAIFJ5X34WLLr2lqlRy1kWwDUaLgCojO79pDYdAqX33zdkPc9hIACILHbPLilrQaD2nXeBwzSAezRcAFAJxhj5rvpJcGDjWtmF37gLBAARxn7xqWSP/CKqXkOpfRe3gQDHaLgAoJJMq7ZS9/6B2n74L9niYoeJACAy2OIi2TKXE5oBF8n4+HET8Y0/AQBwCnyXXxfcDHnXdtmvPjv+CwAgDthFs6WDB/xFYpLMOUPcBgIiAA0XAJwC07CpzDlDA7Wd8I5sfp7DRADgnp31SeDY9DpPJr2GwzRAZKDhAoBTZC7+sZSc7C9y9stO+9htIABwyG7dKK1ZEajNwIvchQEiCA0XAJwic0YdmSEXB2r72QeyuQcdJgIAd+zM4NktNWsttcx0FwaIIDRcAHAazIVXSNXT/EXeYdnPPnIbCAAcsPmHZefMCNRm4EUyxrgLBEQQGi4AOA2merrM8MsDtZ0+QfZgjsNEAFD17DefS6X3saamyfQe4DYQEEFouADgNJkho6S0DH9RkCf72QduAwFAFbKeJzt9YqA25wyVqVbNYSIgstBwAcBpMinVZS4YHajt5xNlc/Y5TAQAVWjFImnXNv+xMTKDR7rNA0QYGi4ACAEzaISUUdNfFBbIThnvNhAAVBFv2oRg0aW3TL2G7sIAEYiGCwBCwKSk+hfQOMLO/ER2/16HiQAg/Oz2zdLKxYHaN/QSh2mAyETDBQAhYgZcJNU8w18UFcpO+a/bQAAQZnZ6mbNbTVtImR2dZQEiFQ0XAISIqVZN5qIrA7WdNUV2b7bDRAAQPjb3kOzsMkvBD7mYpeCBCiS6DuBKQUGBli5dqoULF2rdunXavXu3PM9Tw4YN1adPH40aNUopKSkVvnbWrFmaMmWKtmzZosTERGVmZmr06NFq27ZtFX8WACKNOf8C/5mt/Xul4iLZT96XufZW17EAIOTsV59JhQX+Ir2GTB+WggcqErdnuL766is98cQTmjFjhqy16tKli9q1a6ddu3Zp3Lhx+t3vfqcDBw4c9bo33nhDzz//vDZv3qxOnTqpTZs2ysrK0gMPPKB58+Y5+EwARBKTlCwzYkygtl9+Jrtnl8NEABB6trhItsxiGeb8C2WSkh0mAiJX3J7hSkxM1PDhwzVy5Eg1atQoML5v3z499thj2rBhg/75z3/qrrvuCjy2fPlyTZo0SRkZGXr44YcDr1uzZo0efPBBvfDCC+rQoYPS09Or/PMBEDnMucNkp7wv7c2WSoplJ42TueGXrmMBQMjYeV9I+/f4i8RE/0qtACoUt2e4BgwYoJtvvrlcsyVJZ5xxhn72s59JkubNm6fi4uLAYxMm+H+TM3r06HKvy8zM1LBhw3T48GHNmDFDAOKbSUqSGVnmLNc302V373CYCABCx1or+2lwg3fTd5BMrdoOEwGRLW4bruNp3ry5JKmoqEgHDx6UJBUWFmr58uWSpL59+x71mtKxhQsXVlFKAJHM9B8i1anvL0pKZCe96zYQAITK8oXStk2B0gy/3GEYIPLRcFVg586dkqSEhITA5YHbtm1TUVGRatSooTp16hz1mpYtW0qSNm7cWHVBAUQsk5gkc/GPArWdM5N7uQDEBK/sxu5d+8g0auouDBAFaLgqMHnyZElS165dlZSUJEnKzvYv7VxRsyVJKSkpSktLU25urvLy8qomKICIZvoMLH+Wq+wPKQAQheyGNdKa5YHadwFnt4ATidtFM45l0aJFmjFjhhISEnT11VcHxvPz8yVJycnHXoGnWrVqys3NVX5+vlJTU0/4scaOHXvUWHJysh577DFJUt26dSsbP2QSE/3fGvXq1XOWIdYwp6EXDXN6+KoblfPS/0qS7NfTVPuG25RQ292f7ROJhjmNNsxp6DGnoXeyc7rvH0/pyELwSmrXSXX6sRT8sfB9GnrROqec4Spjy5YtevbZZ2Wt1fXXX68WLVoEHrPWStJxN/QrfQ4AlEodPFK+M440WEWFyv34P24DAcApKt7yvQpmzwzUaZdd6y4MEEU4w3XEnj179Mgjjyg3N1ejRo3SiBHllzctPWNVUFBQ0csl+RfWkHTMDZN/6Mknnzzu49nZ2c6auNLfHOzevdvJx49FzGnoRcuc2qGXSO/9XZJ0+JPxyh8wQia9huNUFYuWOY0mzGnoMaehdzJz6v37Fan055KGTZXTsp0O8jU4Jr5PQ8/VnBpjjlrZvDI4wyUpJydHDz/8sLKzszVw4EBdf/31Rz2n9PK+PXv2VPge+fn5ys3NVVpa2kldTgggfpgBF0rpGf6iIF92+oTjvwAAIozdtU123qxAbUZeJePjx0jgZMT9n5S8vDw9+uij2rp1q3r37q1bb721wssGGzdurKSkJOXk5FTYdG3YsEGS1KxZs7BnBhBdTLUUmaGXBmr7+UTZw7kOEwFA5djJ70ue5y/qN5Lpdb7bQEAUieuGq6ioSH/961+1bt06denSRXfffbd8x/htTXJysjp27ChJmjNnzlGPl4716NEjfIEBRC0zaISUWt1fHM6VnTnZbSAAOEk2e6fsnBmB2oy4SiYhwWEiILrEbcPleZ6eeeYZrVixQu3bt9c999wTWPnkWEaOHClJGj9+vLZv3x4YX7NmjaZNm6bU1FQNHjw4rLkBRCdTPV1m0MhAbad9LHuce0IBIFLYKf+VSkr8RZ36/i0vAJy0uF00Y8qUKZo3b54kKSMjQ6+99lqFz7v++utVo4b/5vbOnTtrxIgRmjx5su6991516tRJJSUlysrKkud5uvPOOwMbJQPAD5mhl8hO+1gqLJAOHpD98lOZoZe4jgUAx2T3Zst+PS1QmxFXypzgF9QAyovbPzGHDh0KHJc2XhW56qqrAg2XJN10001q0aKFpkyZomXLlikhIUEdO3bUFVdcoXbt2oU1M4DoZjJqypx/oey0jyRJ9tPxsgMukjmywToARBo7eZxUXOwvzqgr02+I20BAFIrbhmvMmDEaM2bMKb124MCBGjhwYGgDAYgL5oLLZGdO8v8As3+v7OzpMudf6DoWABzF7tou+9XUQG0uupJfEAGnIG7v4QIAF0ytOjLnDA3Udsp42dJ7IwAggtiP3w7eu1W3gcx5w9wGAqIUDRcAVDFzwWipdEXU3Ttk533hNhAA/IDd8n25v5vMJdfIJHJ2CzgVNFwAUMVMvYYyfQYEajvlv7Kl+9sAQATwPvy3ZK2/aNxMpg/7bgGnioYLABwwF10plW6yvm2TlDXfbSAAOMKu+1ZaGlxQzHfZdTI+9t0CThUNFwA4YBqdKXXtE6i9T96XLf1tMgA4Yq2VN/6N4EDLzHJ/VwGoPBouAHDEd9GVwWL9amnNCndhAECSFs8u93eR7/LrZUrPxgM4JTRcAOCIaZkptescqL0p7ztMAyDe2aJCee/9IzjQqadM+y7uAgExgoYLABzyXXRFsFi+SHbTOndhAMS13AnjpOyd/sLnk++qn7oNBMQIGi4AcKl9V6l5m0Bpp4x3lwVA3CrZv1e5Zc5umUEjZRo1dZgIiB00XADgkDGm3L1cdsHXsru2OUwEIB4deusV2bzD/qJ6uszFP3IbCIghNFwA4Fq3PlKDJv5j63GWC0CVshvXKm/ax4HaXHKNTFqGw0RAbKHhAgDHjC9B5sLRgdrO/lx2/x6HiQDEC1tSIu/N54ObHDdsKjPgQrehgBhDwwUAEcD0HSjVquMviotlp3583OcDQCjYGZOkMov1+K69VSYx0WEiIPbQcAFABDCJSTLDLwvUdtYU2dxD7gIBiHl2727ZD98K1CmDLpIps1UFgNCg4QKACGHOGy6V3jdRkOf/zTMAhIn3n1elgjxJksmooRo33ek4ERCbaLgAIEKYlFSZwSMDtZ0+QbagwGEiALHKLpkjLZkTqDNuulO+mmc4TATELhouAIggZvAoKbmavziUI/vVVLeBAMQce/CAf6GMUplnK7XML3sAhBYNFwBEEJNeQ+b8CwK1/ewD2eJih4kAxBJrrbx/PS8dPOAfSEqW77o7ZIxxGwyIYTRcABBhzLBLpYQjq4Tt3S077wu3gQDEDDv7c2lx8FJCM/oGmUZNHSYCYh8NFwBEGFO7nkzfAYHaTvmvrOc5TAScOmut7OFDstk7ZXdskd3yvezmDf76cK5s6f5PCDu7Z5fsf14JDrTv4r+MGUBYsdECAEQgc8EVst987t+MdPtmKWue1LWv61jAcdm8w9KG1bJrv5XdtE7K3int3S3lHT72i4xPql1XatBEpmETqXEzmdbt/P/38XvhULElJfL+/rSU71+VUKlp8t30K+YYqAI0XAAQgUyjplK3vtKi2ZIkb/L78nXpw30WiDh21zbZhbNlF30jbVzr/yVBpd7Ak/bskvbskl252D8kSalpUut2Mmd3lenSR6Zew5Bnjyf2o7ekNcsDtbnmFpna9RwmAuIHDRcARCjfhVfKO9JwacMa/w9LbTu5DQVIsvl5srNnyH75qbR5w8m/MCFRSkzyHx/Z/+mY8nKl5Qtlly+Uffd1/xmvrn1l+g6U6tEoVIZdMlf2k/cDtel1nkyfge4CAXGGhgsAIpRpeZbUvou0aqkk/1muBBouOGT37Jb9fILsl1P9DVFFqqdJrdrJtGrrP1Nbp4FUt76UllHu8jVbUiIdzpVyc6TdO2R3bJV2bJH9fq2/ibM/uG9x2ybZbZtkJ4/TnrPaK2XARbIde8qkpYfxM45+dtd2/6WEpRqdKXPDLzlbDlQhGi4AiGC+i66Ud6Th0srFshvXyTRv7TYU4o49mCM7eZzszMlSRdsUNGku072/TPe+UuPmJ3VfkElIkDJq+P9r2FSmU8/gx8s/LK1fI/vtUtml86Vtm8q9tui7VSr6bpWUnCzTe4DMoBEyzfhz8UO2oEDei48Fm+NqqfLddp9MSqrbYECcoeECgEjWrrPUvI3/3hhJ9pP3ZW79reNQiBe2qFD20w9kPx0fXGyhVEqqzDlDZQZcKNPozJB+XJNSXerQVaZDV2n0jbK7tssumSM7Z2b5SxgLC2W/murfILxNB/kuukLq1JOzN5KsVyLvtf+TtgTny9x4Z8i/VgBOjIYLACKYMUa+EVf6f0styS76RnbnNpkGjR0nQ6yzq5fLe/M5ade28g+cUVfmgstl+g+RSa1eJVlM/UYywy+Xhl8uu+V7pSydo7zPJ8vm7A8+ae1Kec+ulJq2lLnoCpme58j4EqokX6Sx1sq+85q0pMx+W0Mvka/XuQ5TAfGLhgsAIl3XvlLDJtKOrZK1sp+Ol7nhl65TIUbZw4dk3/+n7JeflX+gerrMiKv8l+8lV3MTTpJp2kI1uvVSxjW3aNeUj2RnTPIvKlNqywbZV5+Q/fDfMhdeIdNvsExSkrO8LtjPPvTPS6nOvWSu/Im7QECcY/MFAIhwxueTufCKQG2/+Vx23x6HiRCr7NpV8h66q3yzZXz+syOPviLfBZc7bbbKMsnV5Os3SAn3PyHffX+VuvQu/4TdO2T/9by8+38ub9pHsoUFboJWMW/OTNn3/xEcaHGWfLf8xn/PHAAnaLgAIAqYPgOkM+r6i5Ji2WkfuQ2EmGI9T94n78v739/5Nyou1bSlfL/7X/muvlmmeuSuBmhat1PCL/8g3wN/k+k9wL+Zcqn9e2XffV3e/bfEfOPlzZkhW3ZFwnoN5bvz/8lUS3GWCQANFwBEBZOYJDP80kBtZ30qm3vQYSLECpt7UN7fHpId/6bkHVmK3eeTuew6+X7/f/7tCaKEadpCvp//j3wPvyhz/oVSYpk7Jw7sK9N4fRxzjZf39XR/s1W6nH56Dfl+9YBMjVouYwEQDRcARA1z3gVSWoa/KMgrf48GcArsji3yHrlHWrE4OFi7nnz3PibfyDEyidF5q7ep30i+62+X79HXZIZeIiUlBx88sE/23dfk3f8LedMnxETj5X0xRfafz0jW+gcyasp3z19kGjZxGwyAJBouAIgaplqKzOBRgdpOnyBbkO8wEaKZXblY3iO/kXZtDw526yvfH5+Rad3OXbAQMrVqy3f1zfI98koFjdde2XdeDTZeRYXugp4i65XI++8bsv96IThY8wz5fvOITJPm7oIBKIeGCwCiiBk8Uiq9H+PQQf/+Q0AleV9MkffMQ8ENcSX/JYS3/U4mLXLv1TpV5RqvIRdLiWVWLQw0XrfImz4xahovm3dY3vOPyE75b3CwVh357nmEvbaACEPDBQBRxKTX8F9aeIT97APZ4iKHiRBNrLXyJo3znxEpvV8ruZp8t93nv4QwxjcMNrVqy/ejn8v3aAWN1/69su+84m+8Po/sxstu3yzvsXulrPnBwSbN5bv3US4jBCIQDRcARBkz7FIp4ci9NXuzZed94TYQooK1Vva9v8t++O/gYK068v32MZnu/d0Fc8DUqnP8xus/RxqvT8fLHs499htVMeuVyPv0A3l/ulvatin4QJfe8t33uEy9hs6yATg2Gi4AiDKmdl2ZvgMDtZ0yXrb0bAVQAeuVyL7xN9mpZbYTqN/Y/0N6s9bugjkWaLweecV/f+QPG6/3/ynv3p/Ke/c12d073AWVZLdvkffX3/n32CpzVttcdKV8t98vk1LdYToAxxOdyw8BQJwzF46W/Wa6f1Wy7ZulpfOkbn1dx0IEsl6J7D//Jjt7RnCwWSv57nqQJcOPMGfUkfnxLbIXXiH7yfuyX34qFRf7HyzIk532sey0j6UO3eQbcIHUuXeVreBo92bLTnxH9utpwctAJf/iGNffIfPDDZ8BRBwaLgCIQqZhU6lbP2nRN5Ik75P35evaJ+bvwUHlWM+TfeO58s1W5tny3fEHmepp7oJFKHNGHZlrfiE74krZzyfJzpoiHT4UfMLKxfJWLpYyasr0PFem93lSq3YyvtBfMGS3b5Gd9Yk/ww/u0zS9B8hcc4tM6TYRACIaDRcARCnfiCvlHWm4tGGNtHqZ1K6z21CIGNbzZN98zn8mtFTbTvLd+f9kSle6RIVMrToyo2+QHXGV7DfT/Xve7dgafMLBA7IzJvnHa9eV6dhDpkM3qV3n01rl0ebsk81a4D+btXbV0U+o11C+q34qw9lsIKrQcAFAlDLN20gdukorl0jyn+VKoOGCjiyQ8c4r/h/cS2WeTbNVSSYlVWbwKNlBI6U1y2VnTZFdNFsqKQ4+aW+27Befyn7xqWR8UpNmMme2lJq2lGncTKp5hlSzlpRew/+450leiXRgn7R7h2z2Tmnzetlvl/kvD65Irdoyo34kc87QqN2MGohn/KkFgCjmu/AKeUcaLq1cIrtxrb8RQ1yzE/4jO2NycKBNB/nu/CPN1ikyxkhtO8m07SSbe1B20WzZ+V9K3y6TbJn7qqwnbfledsv3kmbInu4HbtpS5vzhMv2HylSrdrrvBsARGi4AiGbtOkstM/2XFOrIWa5b73McCi550yfKTngnONAyU767/iiTkuouVAwxaRky5w2XzhvuvwRw+SL/LztWLpEOHjj9D1C7rkznXjLnDpOatea+TCAG0HABQBQzxvjPcr34qH9g0WzZHVv8i2og7nhzZ8m+80pwoNGZ/jNbLBkeFqbGGTL9h0j9h/i3Zti+WXbzemnzBtlN66XsnVLOPqnwGJso16ot1W3g3z+rTQeZdp2keo1osoAYQ8MFANGuax+pYVNpxxbJWtlPP5C58U7XqVDF7LdZsv94JjhQu558dz8kk1HDXag4Ynw+qUlzmSbNpb6DAuPWWqkgTzqYIxkj+RIkn0+qniaTzGWCQDxg42MAiHLG55O58IpAbWfPkN2b7TARqprdtkneC48GF3NIryHfrx+SqV3XbTDIGCOTUl2mXkOZug38G5fXqk2zBcQRGi4AiAGmz/lS6Q/XJcWy0z5yGwhVxh7YJ+9vf5Lycv0DScny/fIPXFYKABGChgsAYoBJTJIZdlmgtl98Knsox10gVAlbUCDvuYelPbv8A8bId/NYmdbt3AYDAATQcAFAjDDnDZfSM/xFQX75ZcERc6y1sm/8Tfr+u8CYueqnMt37O0wFAPghGi4AiBGmWorMkIsDtf18gmxBvsNECCc7Zbx/L6gjzKARMkMvcZgIAFARGi4AiCFm0Eip2pH9lg4dlP3yM7eBEBZ22ULZD94MDrTvInP1z1lOHAAiEA0XAMQQk5YhM+CCQG0/+1C2uMhhIoSa3bFV3qtPSNb6B+o2kO+W38gkJLgNBgCoEA0XAMQYM/RSKeHINov7smW/+dxtIISMzTss7/m/BFckrJYi3x2/l0lnry0AiFQ0XAAQY8wZdWT6Dw7UdtI42SLOckU763nyXn/Sv8H1Eb6f3C3TtIW7UACAE6LhAoAYZEaOCZ7l2rtb9uupbgPhtNmP35aWzgvUZuQYmR6sSAgAkY6GCwBikKlTX+a8YYHaTnpPtqjQYSKcDrvwa9lJ44IDXXrLXHKNu0AAgJNGwwUAMcpcdJWUeOQs1/49sl986jYQTondtkneP54JDjQ6U76fjZXx8U84AEQD/rYGgBhlateVOf/CQG0/eV+2oMBhIlSWLciX99LjUul+atXT/ItkpFZ3GwwAcNJouAAghpmLrpSSkv3FgX2ysya7DYSTZq2V/feL0vbNgTHfT8fKNGjsMBUAoLJouAAghplatWUGXhSo7ZTxsvl57gLhpNmvpsrOmRGozYVXyHTp5TARAOBU0HABQIwzF14hJVfzFwcPyM6Y5DYQTshuWi/79svBgbM6yFx2nbtAAIBTRsMFADHO1KglM3hUoLaffiCbd9hhIhyPzTss7+XHpeIje6dl1JTv57+RSUhwGwwAcEpouAAgDpjhl0vVUv1F7kHZ6R+7DYQKWWtl33hW2rXdP2CMfDePlTmjjttgAIBTRsMFAHHAZNSQGXpxoLaffSR7+JDDRKiInTFJduHXgdqMulqmQzeHiQAAp4uGCwDihBl2mVS6nHheruxUznJFErvhO9lxfw8OtO8iM+pqd4EAACFBwwUAccKkpcsMvTRQ22kfyR7McZgIpWzuIf99WyXF/oGatf2XEvq4bwsAoh0NFwDEETP0Eql6ur/Iz5OdPM5tIMhaK+8fT0t7dvkHjE++W+6RqXGG01wAgNCg4QKAOGKqp8mMuDJQ2xmTZXfvcJgI9rMPpaXzArW5/DqZzI7uAgEAQoqGCwDijBk8Sqpd11+UFMt++JbbQHHMrl0pO/6N4ECnnjIXjHYXCAAQcjRcABBnTFKyzKXXBmo7b5bsxnUOE8Une/CAvJf/V/I8/0DtuvL99G4ZH/80A0As4W91AIhDpu9AqUnzQO2VPcuCsLOeJ++1J6X9e/wDCQny3XKvTHoNt8EAACFHwwUAccj4EuS74sbgwMolsisWuwsUZ+zk96SVwfk2V94k07qdw0QAgHBJdB0gGhUWFurDDz/U119/rezsbKWnp6tLly66+uqrVadOHdfxAODkdOwhte0krV4mSfLe+7vseYNlEvinIZzsqqWyH/8nONCtr8yQS9wFAgCEFWe4KqmwsFB//vOf9f777ys/P189e/ZUnTp1NHPmTP32t7/Vjh2s9gUgOhhj5LvqJ5Ix/oGtG5XHZshhZffvlffa/0n2yH1b9RrKd9OvZEq/BgCAmEPDVUkffPCBVq9erczMTD3zzDP69a9/rUceeUQ33HCDcnJy9OKLL7qOCAAnzTRvI9NvcKA++Par8g4ddJgodtmSEnmvPiHl7PcPJCbJ94vfypTuiwYAiEk0XJVQXFysKVOmSJJ+9rOfKSUlJfDYqFGj1Lx5c61atUrr1693FREAKs1cfp1Uzf/3mc3Zr0Pv/cNxothkP/6PtGZ5oDY/+rlM89YOEwEAqgINVyV8++23ys3NVYMGDdSyZcujHu/Tp48kacGCBVUdDQBOmalVR+ai4GbIhye9J7tzm8NEsadg4WzZyeMCtek9QOb8CxwmAgBUFRquSti4caMkVdhsSVKrVq3KPQ8AooUZdqlUu56/KC6W997f3QaKISW7d2r/0w8FBxo2lbn+du7bAoA4QcNVCdnZ2ZJ0zJUIa9euXe55ABAtTHI1mStvCg4snSebNd9Znlhhi4u0/4k/yB484B9ITpbv1t/KpKS6DQYAqDKs/VsJ+fn5kqRq1apV+HjpPV2lzzuRsWPHHjWWnJysxx57TJJUt27dU4kZEomJ/m+NevXqOcsQa5jT0GNOQ8tedLn2fTNNhcv9+0OZca+r7rmDZaqlnOCVOJacV/5Ph1cH79uqeeu9Su3a02Gi2MCf/dBjTkOPOQ29aJ1TznBVgrX2tB4HgEhmjNEZt94r+RIkSSU7t+nQf990nCp65c2cosOT3w/UqUNGKXXwSIeJAAAucIarElJT/ZeAFBQUVPh46XjZ1QuP58knnzzu49nZ2c6auNLfHOzevdvJx49FzGnoMaehV69pC1W/5God/vBtSVLu+H8pr3MfmQaNHSeLLnbTenkvPBqoE1u1VcHoG/leDRH+7Icecxp6zGnouZpTY4waNWp0yq/nDFcllF7it2fPngof37t3b7nnAUA0Sr/6Z1KtI/eqFhfL+8/LnMGvBJt7UN6Lj0qFhZIkk1FDtX77qExyxZejAwBiGw1XJTRv3lyStGHDhgofL91/q/R5ABCNfKnV5bv6Z8GBFYtl53/pLlAUsZ4n77Unpeyd/gFjVGvsn5TY4NR/MwoAiG40XJXQrl07Va9eXTt37qyw6Zo7d64kqXv37lUdDQBCq8c5UodugdL+55XgSns4JjvhHWn5wkBtLr1W1br1cZgIAOAaDVclJCYm6sILL5Qk/f3vfy+3GuHEiRO1ceNGtWvXTm3atHEVEQBCwhgj33W3SaWXwR3KkX3nNbehIpxdOl924jvBga59ym0oDQCITyyaUUmjR4/WsmXLtHr1at11111q166dsrOz9d133ykjI0O3336764gAEBKmXkOZy6+Tffd1SZKdN0u29/kyXXo5ThZ57K5t8l4vsxBSgyby/eRuGR+/1wSAeMe/BJWUnJysBx54QFdccYWSk5M1f/587dq1SwMGDNDjjz+uhg0buo4IACFjBo+SWrcL1N6/X5A9nOswUeSx+Yf9KxLmHZmXainy3fY7meppboMBACICZ7hOQXJysq6++mpdffXVrqMAQFgZX4J8N94p7093ScXF0v49suNek7npLtfRIoL1SvyLZGzdGBgzN/5Kpkkzh6kAAJGEM1wAgOMyjc6UGfWjQG2/ni678BuHiSKHHf8vaem8QG0uGC1fr3MdJgIARBoaLgDACZkLRkstMwO19+Zzsvsq3pMwXnjfTJf9dHxwoEtvmdHXuwsEAIhINFwAgBMyiYny3TxWqpbiHzh8SN4/npb1PLfBHLFrVsj+6/ngQJPm8t08VsaX4C4UACAi0XABAE6Kqd9Y5kc/Dw6sWio77SN3gRyx2zfLe/5h/z1tkpRRU747/59MSnW3wQAAEYmGCwBw0sw5Q6Xu/QO1Hf8v2XXfOkxUtez+vfKeeUgqXakxMUm+2++XqVPfbTAAQMSi4QIAnDRjjHw33CHVquMfKCmW99Ljsjn7neaqCjb/sLxn/yzt2eUfMEa+m/9Hpk17t8EAABGNhgsAUCkmLUO+W34jJRy5X2n/Hnmv/K9sSYnbYGFkiwr9e21tWhcYM2N+JtOj/3FeBQAADRcA4BSYszrIXPXT4MDqZbIfvOkuUBjZkhJ5rzwhrVoaGDNDL5Vv6CUOUwEAogUNFwDglJjBo2R6DwjU9tMP5M2d5TBR6FnPk/3nM9KSOYEx0/t8mat+4jAVACCa0HABAE6JMUbmhjukJs0DY/afz8iuXu4wVehYa2Xffkl2zszgYJfeMj+5W8bHP58AgJPDvxgAgFNmqqXId/vvpPQM/0BxsbwX/iK7bZPbYKfJep7sv1+QnTUlONius3y/uFcmMdFdMABA1KHhAgCcFlO/sXy//H9SUrJ/4HCuvGcekt2/122wU2S9Etl//k32i0+Dgy0z5bvj9zKlnyMAACeJhgsAcNpM63by3TxWMsY/sHe3vKcfkD14wG2wSrIlJbKvPyU7+/PgYOt28t39kExKqrtgAICoRcMFAAgJ072/zJifBQe2bpT3f3+ImqbL5h+W99yfZed9ERzM7Cjf3Q/KVE9zFwwAENVouAAAIeMbeonMiKuCA6VNV4RvjGz37ZH3199JyxcFBzt0le9XD8ikVHcXDAAQ9Wi4AAAhZS67TmbEmODA1o3ynvi97J7d7kIdh92yQd6jv5E2bwgOdu8v3y//IFOtmrtgAICYQMMFAAgpY4zMZdfKjCzTdG3fLO+R/5HdsMZdsAp430z3N1v7sgNjZvhl/tUIWSADABACNFwAgJAzxshceq3MxT8KDubsl/e/98ub/5W7YEfYwgJ5bzwr+49npMJC/6DxyVzzC/mu+in7bAEAQobNRAAAYWGMkbnkGnln1JV960WppEQqKpR95a/yvl/jv/TQwVkku3GdvH8+I235PjiYliHfz8bKdOpR5XkAALGNhgsAEFa+84bL1m8k74VHpcOHJEn2sw9lVyyW76e/lmnWqkpy2IIC2Qlvy079SPK84AMtM+X7xW9l6tSrkhwAgPjCNRMAgLAzbTvJd/8TUsOmwcGtG+U9co+8j9+Wzc8L28e2nie78Bt5D90p++kH5ZotM+Ri+e59lGYLABA2nOECAFQJ06CxfH94SvaDN2WnT/APlhTLTnhHdtYUmVFXy5w3XCYxKSQfz1orZS2Q99G/y69AKEn1Gsp3/R0y7buE5GMBAHAsNFwAgCpjqlWT+dHPZTv3kvfPvwVXB8zZL/v2y7KffiBz7jCZ/oNlap/aWSebs1927izZr6dJWzeWf9Dnkxl2mczFP2bJdwBAlaDhAgBUOdOhq3wPPis75b+y0z8OrhS4Z5fsR2/Jfvy21L6LTIduMq3bSs3bHHOBDVtUJG1cK7t2pezq5dKqJf4FOn6oW1/5Lr1WpknzsH1eAAD8EA0XAMAJUz1NZvQNsoNH+i8r/Gpq8P4qa6WVS2RXLpGVpIREqXZdqXq6VD1N8vmkgznSoQPSgf1SSfGxP1Cnnv5Gq3nrKvisAAAoj4YLAOCUqVVH5vo7ZEeOkZ09Q/ab6dKu7eWfVFIs7d5x8m96Rl2ZfoNk+g2WadgktIEBAKgEGi4AQEQwtevJjBwjO+Iq6buVsssWyK5fLX3/nVRYcPwXJydLLTJl2nSQaddJattRxpdQNcEBADgOGi4AQEQxxkiZZ8tkni1JssXF0o7N/oU1cnP9e3mVFEsZNWUyakoZtaQGjUK2uiEAAKFEwwUAiGgmMVFq2tJ/7DgLAACVxcbHAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECaJrgPg2IwxriNERIZYw5yGHnMaesxp6DGnocechh5zGnrMaehV9Zye7scz1loboiwAAAAAgDK4pBAAAAAAwoSGCxW67777dN9997mOEVOY09BjTkOPOQ095jT0mNPQY05DjzkNvWidU+7hQoUKCwtdR4g5zGnoMaehx5yGHnMaesxp6DGnocechl60zilnuAAAAAAgTGi4AAAAACBMaLgAAAAAIExouAAAAAAgTNiHCwAAAADChDNcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmNFwAAAAAECY0XAAAAAAQJjRcAAAAABAmia4DoGrk5+dr3rx5Wrt2rb777jtt3LhRxcXFuuaaa3TZZZcd97V79uzRu+++q6VLl+rQoUOqW7eu+vfvr8svv1zJycmVzrJlyxaNGzdOK1asUH5+vho2bKhBgwZpxIgR8vmi/3cADz74oFauXHnc5xhj9O67757U+82cOVMvvPDCMR/v37+/7r777spEjDorVqzQQw89dMzHzzrrLP3lL3+p9PsuXLhQH3/8sb7//ntJUosWLXTJJZeoR48epxo1amzdulXz589XVlaWtm/frgMHDigtLU1t27bVyJEj1b59+0q9Xzx9nxYWFurDDz/U119/rezsbKWnp6tLly66+uqrVadOnUq9V25urt577z3NmzdP+/fvV61atdSrVy+NGTNGaWlpYfoMIkdBQYGWLl2qhQsXat26ddq9e7c8z1PDhg3Vp08fjRo1SikpKSf9fnfccYd27959zMefeuopNWnSJBTRI9qJ/h26//771bVr15N+v3j/Pj3Rv0GlxowZoyuvvPKEz4un79P169crKysr8PPnvn37lJSUpLfeeuu4r5s1a5amTJmiLVu2KDExUZmZmRo9erTatm1b6Qye5+mTTz7R559/rh07diglJUVnn322xowZo6ZNm57qp1YpNFxxYseOHXruuedO6XV/+MMflJOTozPPPFPt2rXT+vXr9d///lfLli3TAw88oKSkpJN+vzVr1ujPf/6zCgoK1KZNG9WrV0+rVq3Sm2++qdWrV2vs2LEyxlQ6ZyTp2rWr6tWrV+Fj69ev1+bNm9WuXbtKv2/z5s3VokWLo8bPOuusSr9XtGrQoEGFc9egQYNKv9fkyZP1z3/+UwkJCerUqZMSExOVlZWlxx9/XDfddJNGjBgRisgR689//rP27t2r1NRUnXXWWcrMzNSWLVs0b948zZ8/XzfccINGjhxZ6feN9e/TwsJC/fnPf9bq1at1xhlnqGfPntq9e7dmzpypRYsW6eGHH1bDhg1P6r0OHjyoP/zhD9q+fbsaNGigXr16acuWLfrkk0+0ePFi/eUvf1FGRkaYPyO3vvrqK7388suSpDPPPFNdunRRXl6e1qxZo3Hjxunrr7/Wgw8+qJo1a1bqfQcMGFDhePXq1U87czTp06dPhQ1r7dq1T/o9+D6VatWqdczvKc/z9OWXX0pSpf9tj4fv0/fff18LFiyo1GveeOMNTZo0ScnJyercubOKioqUlZWlpUuXauzYserdu/dJv5e1Vk8//bTmzJmjtLQ0de/eXQcPHtTcuXO1aNEiPfDAA1Xy7xMNV5xISUnR4MGD1aZNG7Vu3Vpz587V+PHjT/i6F198UTk5Obrooov0k5/8RJJUUlKip556SvPmzdMHH3ygMWPGnFSGkpISPfvssyooKNANN9ygUaNGSfKffXv44Yc19/+3d+9BUZVvHMC/y4IiCMKAIAgpYoiEgqGYhaWmjhfSwttUUqZOOV7G7d6kjVn+UnNyspSaNFLxMtnFUiMammiSCjEQg0RuAoHgIgRCC8iy7O8P5pxc2GV3YQ+wy/fzl5zLy7vHh3ff55z3vO+FC/j5558xc+bM7n/QfqCrJ4avv/46AODBBx80u1zhbuJAFhwcjA0bNvS4nIqKCiQkJMDBwQHbtm1DUFCQuP2NN95AQkICJk2aBB8fnx7/rv7Kz88PK1euxH333Qd7+/++CpKTk3Hw4EEkJCQgLCzM7Lt/th6np0+fRl5eHoKCgrB161axM3vu3DkcPXoUH330kUl3woH2TkVlZSUiIyPx/PPPQy6XAwDi4+ORlJSEI0eOYOPGjZJ9lv7A3t4ec+fOxcKFC3X+3mpra7Fr1y4UFxfj8OHD2Lx5s1nlWqKdsAWxsbHw8vLqURmMU2DkyJEGY+rSpUs4f/48PDw8EBISYla5AyFOg4KCMHr0aAQGBiIwMBDPPvtsl8fn5OTgu+++g4uLC3bs2CG2C/n5+XjzzTcRFxeHkJAQDB061KTfn5KSgrS0NPj4+GD79u1wc3MDAKSlpWHv3r344IMP8P7774txLRXrH79FJhkxYgTWrVuH2bNnIyAgwKShe4WFhcjNzcWwYcOwcuVKcbtcLsfatWshl8vx/fffo7W11aQ6pKenQ6lUYtSoUWKyBbQng2vWrAHQ3mmxVZWVlSgsLISDgwOmTZvW19UZ0BITE6HRaDBnzhwx2QIAX19fPPbYY9BoNPj+++/7sIbS27p1K6KionSSLQCYM2cOwsLC0NbWht9//72Patc/tba2IikpCQCwZs0anScH0dHRGDVqFHJzc3Ht2jWjZdXV1eH8+fM67akgNjYWrq6uSE1NRV1dncU/R3/y0EMPYe3atZ1ubri7u4vfC+np6SZ/z5BlMU6NE55uTZ8+3SZei7C0Rx99FMuXL0dERISY7HTl7NmzAICYmBiddiEoKAhz5sxBY2MjUlJSTP79Qr/yySef1Pn99913HyZPngylUomLFy+aXF53MTLIoMzMTABAREREp2GDbm5uGD9+PFQqFfLy8kwqLyMjA0B7kHcUEBAAb29vlJWVoaqqqoc175+ERjkiIsKmhgtYIyG29cWikAwL8ToQjRo1CkD7Uwb6z9WrV6FSqeDt7Y2AgIBO+6dOnQoAJg2fuXTpErRaLUJCQjp1QhwcHBAREYG2tjZkZWVZoupWSYhDtVqNhoaGPq7NwMQ47Vpzc7PYWZ8+fXof18b6tbS0ICcnB4D+72dhm6nfz1VVVSgvL8egQYNw77339ri8nuCQQjJImEhAX8dC2J6Tk4PS0lLcc889RssrLS01Wp5SqURpaWmPh0D0R0LC1Z3hhED7+18JCQloamqCm5sbQkNDzR6+YO1u3LiBEydOoKGhAS4uLggODkZ4eLhZdxVVKhWqq6sBQO+7Rh4eHnBxccHNmzfR2Ng4IJNjpVIJACbdjezIluPUWBs2ZswYneN6UlZAQABSUlLEdnggEuJQLpebPHxIcObMGdy4cQMODg7w9/dHZGQkXF1dpahmv/bTTz/h33//hUwmg4+PDyIjI+Hp6Wny+YzTrqWnp+P27dsICAiAv7+/2eczTnVVVFRArVbD1dVV7wREQhya0sYC//Vj/f39O43m6E55PcGEiwyqqakBAIOzbgnbhc6rMcJxhsoTXuI1tTxrkp+fD6VSCRcXF7NmhrpTZmam+GQGaH8RNSQkBAqFolsdY2uUl5fX6YnqXXfdhRdffNHk962E+HJ2djY4+5mHhwcaGhpQXV2Nu+66q2eVtjI3btwQ42zy5Mlmn2/LcWrJNkw4xtDkBea2r7YoMTERQPtEROZMzgQAx44d0/n5yJEjeOaZZzBr1iyL1c8adHxXOyEhAUuWLDFpJj2AcWrMncMJu4NxqstYG+vo6AhnZ2eoVCo0NTVhyJAhPSqvN+OXCRcZ1NzcDAAGp34fPHiwznGmliec15HQ+TW1PGvyyy+/AGifGlvfXZauuLm5YdmyZZgyZQq8vLzQ0tKCwsJCHD9+HFeuXMGuXbvwzjvv2PTYcScnJyxatAhTp04VE6uSkhKcPHkSBQUF2LFjB/bs2WPS0yhjcXjnPluMxa5oNBrExcVBrVbj/vvvF5/YmGIgxKkl2zBjZQnbb9++bXY9bUFmZiZSUlIgl8uxYsUKk8+LiIhAaGgoxowZA1dXVyiVSqSkpCAxMREff/wxhg4datYMZ9Zq/PjxmDVrFsaNGwd3d3dUV1cjLS0NX3/9NU6dOgUnJyeTZmJlnBpWV1eH7Oxs2NnZISoqyqxzGaf6Get3Au0xp1Kp0NzcbDThsnQ/tieYcFmJ9957D2VlZWads3HjRowdO7bbv1Or1QKAwWnahf2WYunyusvS17q1tVWcfKA7wwnDw8N1noo5OTlh8uTJCA0Nxauvvopr167ht99+M7vB7009vaYBAQGdhrSEhobi7bffxvbt25Gbm4ukpCTExMQYLddYXFsLKdqE+Ph4XL16Fd7e3li7dq1ZZdtCnBpjrI0ypw2zlTiUQnl5OT788ENotVrExsbqHfpryOrVq3V+9vf3x1NPPQVfX1988sknOH78+IDoyHZMUn19fRETE4PAwED873//w6lTpzB79myja2kyTg1LTU1FW1sbwsPDzX56zzjVz5R4605fsT/ELxMuK3Hz5k1UVFSYdU5P7zgJdw4MldPS0gIAJi9K6ejoCJVKZbA8Ybs5i1xKwdLXOisrCw0NDfDx8bHoWg+Ojo6YP38+4uPjkZWV1a87slLFr52dHRYvXozc3FxcvnzZpIRLiOuu7mj1l1jsiqWv6Zdffonk5GQMGzYMW7ZsMfudGUOsKU6NMdYmmhM3xuJQKKurJ7G2qKamBu+88w5UKhWio6Mtth7erFmz8Pnnn6OyshJVVVU2+Z6wKcLCwhAYGIiioiLk5+cjNDS0y+MZp4b19L1sfQZ6nBprYwHz+p7CMf2h38mEy0rs2rWr13+nh4cHiouLxXe5OhK2m/oCrqenJ1QqFWpqasTZp+70zz//mFWeVCx9rYXhhFLMYCQMr+vvU/JKGb/CIrOmXgMhvoQhCfoaWnNjuy9Y8pomJSWJw4y2bNli8sK9prKWODVGiAdDbaI5bZhwjHBOR9YQg5ZWX1+PHTt2oLq6GjNmzEBsbKzFyrazs4O3tzdu3bqF2traAdeRvdOIESNQVFRk0t8j41S/8vJyFBcXw9HREVOmTLFYuQM9To21sc3NzVCpVHB2djY6nNCU8nozfq13MD1JThjGUVxcrHe/sN3USQWEJMtYefqSMWvV2NgoTjcqRcKlUqkA9O8nMVIz9xo4OzuLjau+mbVqamrQ0NAAT0/PATFD4fnz5/HZZ59h8ODBeO2118wavmUqW4lTY22YsP6WKW3YQGwPu9LU1ISdO3fi+vXriIyMxLp16yw+DMhW4rCnzLkOjFP9hBupkZGRFn+6N5Dj1NfXFw4ODqivr9ebJJnb7xS+z8rKyvSu5WdueT3BhIsMEtYsyMjIgFqt1tlXV1eH3NxcODk5ITg42Kzy0tLSOu0rLi6GUqmEn5+fTd3RSUtLg1qtxrhx4+Dt7S1J+QDMmtzA1ly4cAGA4WmL9ekqFoX37fSt2WFrMjMzERcXB7lcjpdeesnkv2Vz2UqcBgcHw8nJCUqlUm8HVIhFU2InPDwcMpkMubm5uHXrls4+tVqNjIwMyGQyTJo0yTKV78fUajXeffddFBUVISwsDAqFwuKTq5SVlaGiogKDBw/GyJEjLVq2Namvr0dubi4A09pMxmlnWq0Wv/76KwDLDicEGKeDBg0Sh7nq+34WtkVERJhUnpeXF0aOHImWlhad2XO7W15PMOEig8aOHYtx48bh1q1bOH78uLhdo9Hg0KFD0Gg0mDdvXqdZ9/bv3w+FQoH09HSd7ZGRkfDy8kJpaam48jfQ/oj4008/BQBER0dL+Il6nzljvBUKBRQKRaehG4mJiZ3Gz7e2tuKLL75AWloaBg0ahBkzZliszv1RcnJyp4VPtVotkpOT8d1330Emk2Hu3LmdzjN0TRcsWAA7OzskJycjPz9f3F5ZWYnTp0/Dzs7OYu+O9FdXr17F3r17AbRfp7CwMJPOG8hxam9vj3nz5gFon2Dkzs977tw5lJaWIjg4WGdikqSkJCgUCpw4cUKnLHd3dzzwwANobW0V21PBsWPHUF9fj6ioKKufSt+YtrY27Nu3D3/99RfGjx+Pl156yehMroauaVZWlviU8U6lpaXYu3cvtFotZs2aZfZMsdYmPz8fOTk5nSYXqKqqwp49e3D79m1MnjxZZ6psxqnpcnNzcfPmTbi7u3f5DhzjtHsWLlwIoH1Jg8rKSnF7fn4+fvzxRwwZMqTTtPmFhYVQKBR46623OpUn9CuPHz+uc9PgwoUL+OOPP+Dl5WXRYaGGDMz/zQFqz5494pht4VHtDz/8IK6S7ubmhpdfflnnnPXr12Pr1q1ITExETk4O/Pz8UFRUBKVSibvvvlvvJAXV1dWoqKhAY2OjznZ7e3ts2rQJb7/9No4ePYrff/8dnp6euHr1KmprazFlyhSr75Dd6Z9//sGVK1dgb2+PadOmGT1emACh42Pvw4cP48SJE/Dz84OnpyfUajVKSkpQW1sLBwcHbNq0yeAaKbbim2++QXx8PPz8/DB8+HAAwN9//42qqirIZDKsWrVK79MTQ9fU19cXK1euxNGjR7Ft2zZMnDgRcrkcf/75J1paWsTZomzZ7t270dLSAi8vL1y8eFFsB+4UHByMhx9+WGfbQI/TmJgYZGdnIy8vD5s3b0ZwcDCqq6tRUFAAFxcXrF+/Xuf4+vp6VFRUoLa2tlNZq1atQkFBAS5cuACFQoHAwECUlZWhrKwM3t7eePrpp3vrY/WZpKQk8eaci4sLDh06pPe42NhYcUFYQ9c0Pz8fX375JYYPHw5vb2+4urqiqqoKxcXF0Gg0CAkJwRNPPCHtB+oHKioqEBcXB3d3d/j4+MDNzQ01NTW4du0a1Go1/P398dxzz+mcwzg13Z1rb3X1JJZx2i4zMxNfffWVzrbW1lZs2bJF/HnJkiXiyICJEydiwYIFSExMxCuvvIIJEyZAo9Hgzz//RFtbGzZt2tRpUqfbt2+LiyZ3NHPmTFy6dAnp6elQKBSYMGECGhoacOXKFfG7qTeSWyZcA0hJSQlu3ryps62mpkZMvoSO7J18fHywe/dunDp1CllZWUhPT4eHhwdiYmIQExNjdErZjsaNG4edO3fi1KlTuHLlCkpKSuDt7Y3o6GgsXLjQqtfo6ej8+fPQarW49957ezTj29KlS5Gfn4/r16+jvLwcWq0WHh4emD17NqKjo20+MQDa71BdvnwZ5eXlyM7Ohkajgbu7O6ZPn4758+d3a/mD6OhojBgxAmfPnhWH2IwZMwaLFi3q1oK/1kZ4T6CqqgpVVVUGj+uYcBkyUOJ00KBB2LZtG06fPo3U1FRcvHgRzs7OeOihh7BixQqzXr52dXUV28OLFy8iPT0dw4YNw7x587B8+XKLzRTZn/3777/ivzuOirjTsmXLxITLkPDwcNTU1KCoqAilpaVobGzEkCFDEBwcjKioKMycOdOmvmMMGTt2LObOnYuCggKUl5cjLy8PgwcPxujRozFt2jTMnTvXrO9uxul/1Gq1OAytu+9lD7Q4ra+vR0FBgc42rVars62+vl5n/6pVqzB69GgkJSUhOzsbcrkcoaGhWLJkidlD3+3s7PDCCy8gMTERKSkpyMjIECc7WbFiBfz9/bv/4cwg0/aXxY+IiIiIiIhsjO2k0ERERERERP0MEy4iIiIiIiKJMOEiIiIiIiKSCBMuIiIiIiIiiTDhIiIiIiIikggTLiIiIiIiIokw4SIiIiIiIpIIEy4iIiIiIiKJMOEiIiIiIiKSCBMuIiIiIiIiiTDhIiIiIiIikggTLiIiIiIiIokw4SIiIiIiIpIIEy4iIiIiIiKJMOEiIiIiIiKSCBMuIiIiIiIiiTDhIiIisoBvvvkGy5cvx+OPP47CwkK9x2RmZmLFihVYvnw5UlNTe7mGRETUF5hwERERWcDixYsxYcIEaDQa7Nu3D01NTTr7a2trERcXB61WiwcffBBRUVF9VFMiIupNTLiIiIgsQCaTYdOmTRg2bBiUSiUOHjwo7tNqtdi/fz/q6+sxYsQIrF27tg9rSkREvYkJFxERkYW4ublh/fr1kMlkSE1Nxc8//wwA+Pbbb5GdnQ25XI7NmzfD0dGxbytKRES9hgkXERGRBU2aNAkLFy4EAMTHx+OXX37B559/DgB4/PHHERgY2JfVIyKiXibTarXavq4EERGRLWltbcXWrVtx7do1cVtYWBhef/11yGSyPqwZERH1Nj7hIiIisjB7e3usX79e/NnJyQkbNmxgskVENAAx4SIiIpLAjz/+KP67qakJJSUlfVcZIiLqM0y4iIiILCwjIwNJSUkAgFGjRkGr1eLAgQOoq6vr24oREVGvY8JFRERkQcJ6WwAwY8YMbN++HcOHD8etW7dw4MAB8NVpIqKBhQkXERGRhbS1tWH//v1oaGiAj48PVq9eDScnJ2zevBlyuRyXL1/GuXPn+rqaRETUi5hwERERWciZM2f0rrcVFBSEpUuXAgBOnjypM3shERHZNiZcREREFlBYWKiz3taYMWN09j/22GO455570Nrain379qG5ubkvqklERL2MCRcREVEPNTU1Yd++fdBoNJg4cSIeeeSRTsfY2dlh48aNGDp0KCorKxEfH98HNSUiot7GhY+JiIiIiIgkwidcREREREREEmHCRUREREREJBEmXERERERERBJhwkVERERERCQRJlxEREREREQSYcJFREREREQkESZcREREREREEmHCRUREREREJBEmXERERERERBJhwkVERERERCQRJlxEREREREQSYcJFREREREQkESZcREREREREEmHCRUREREREJBEmXERERERERBJhwkVERERERCQRJlxEREREREQSYcJFREREREQkESZcREREREREEvk/I/QLzUNYSZsAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : x**2 + 10*np.sin(x)\n",
"fp = lambda x: 2*x + 10*np.cos(x)\n",
"\n",
"x = np.arange(-10, 10, 0.1)\n",
"\n",
"plt.plot(x, f(x))\n",
"plt.xlabel('x')\n",
"plt.ylabel('y')\n",
"plt.legend([r\"$x^2 + 10\\sin(x)$\"])"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "757e1afe",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = -2.3787, err = 1.621e+00\n",
"x1 = -2.4832, err = -1.045e-01\n",
"x1 = -2.4795, err = 3.666e-03\n",
"x1 = -2.4795, err = 4.253e-06\n",
"x1 = -2.4795, err = 5.736e-12\n"
]
}
],
"source": [
"# 초기값 -4\n",
"x0 = x = -4\n",
"\n",
"Err = []\n",
"itmax = 5\n",
"for i in range(itmax):\n",
" # Newton-Raphson\n",
" xn = x - f(x) / fp(x)\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "cd14fbd5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = -0.2717, err = -1.272e+00\n",
"x1 = 0.0154, err = 2.872e-01\n",
"x1 = 0.0000, err = -1.541e-02\n",
"x1 = 0.0000, err = -2.251e-05\n",
"x1 = 0.0000, err = -5.067e-11\n"
]
}
],
"source": [
"# 초기값 1\n",
"x0 = x = 1\n",
"\n",
"Err = []\n",
"itmax = 5\n",
"for i in range(itmax):\n",
" # Newton-Raphson\n",
" xn = x - f(x) / fp(x)\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "b55da5b9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 6.5628, err = 4.063e+00\n",
"x1 = 4.5472, err = -2.016e+00\n",
"x1 = 3.0957, err = -1.451e+00\n",
"x1 = 5.7397, err = 2.644e+00\n",
"x1 = 4.3537, err = -1.386e+00\n",
"x1 = 2.5082, err = -1.846e+00\n",
"x1 = 6.5195, err = 4.011e+00\n",
"x1 = 4.5492, err = -1.970e+00\n",
"x1 = 3.1005, err = -1.449e+00\n",
"x1 = 5.7449, err = 2.644e+00\n",
"x1 = 4.3563, err = -1.389e+00\n",
"x1 = 2.5186, err = -1.838e+00\n",
"x1 = 6.4670, err = 3.948e+00\n",
"x1 = 4.5496, err = -1.917e+00\n",
"x1 = 3.1014, err = -1.448e+00\n",
"x1 = 5.7459, err = 2.645e+00\n",
"x1 = 4.3568, err = -1.389e+00\n",
"x1 = 2.5206, err = -1.836e+00\n",
"x1 = 6.4575, err = 3.937e+00\n",
"x1 = 4.5495, err = -1.908e+00\n",
"x1 = 3.1010, err = -1.448e+00\n",
"x1 = 5.7455, err = 2.644e+00\n",
"x1 = 4.3566, err = -1.389e+00\n",
"x1 = 2.5197, err = -1.837e+00\n",
"x1 = 6.4618, err = 3.942e+00\n",
"x1 = 4.5495, err = -1.912e+00\n",
"x1 = 3.1012, err = -1.448e+00\n",
"x1 = 5.7457, err = 2.645e+00\n",
"x1 = 4.3567, err = -1.389e+00\n",
"x1 = 2.5201, err = -1.837e+00\n",
"x1 = 6.4596, err = 3.939e+00\n",
"x1 = 4.5495, err = -1.910e+00\n",
"x1 = 3.1011, err = -1.448e+00\n",
"x1 = 5.7456, err = 2.645e+00\n",
"x1 = 4.3566, err = -1.389e+00\n",
"x1 = 2.5199, err = -1.837e+00\n",
"x1 = 6.4606, err = 3.941e+00\n",
"x1 = 4.5495, err = -1.911e+00\n",
"x1 = 3.1011, err = -1.448e+00\n",
"x1 = 5.7456, err = 2.645e+00\n",
"x1 = 4.3567, err = -1.389e+00\n",
"x1 = 2.5200, err = -1.837e+00\n",
"x1 = 6.4601, err = 3.940e+00\n",
"x1 = 4.5495, err = -1.911e+00\n",
"x1 = 3.1011, err = -1.448e+00\n",
"x1 = 5.7456, err = 2.645e+00\n",
"x1 = 4.3566, err = -1.389e+00\n",
"x1 = 2.5200, err = -1.837e+00\n",
"x1 = 6.4604, err = 3.940e+00\n",
"x1 = 4.5495, err = -1.911e+00\n"
]
}
],
"source": [
"# 초기값 2.5\n",
"x0 = x = 2.5\n",
"\n",
"Err = []\n",
"itmax = 50\n",
"for i in range(itmax):\n",
" # Newton-Raphson\n",
" xn = x - f(x) / fp(x)\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "markdown",
"id": "28a1de15",
"metadata": {},
"source": [
"### 문제점\n",
"- 수렴하지 않을 수 있다.\n",
" * Local extrema 근처에서 진동할 수 있다.\n",
"\n",
"- 초기 추정 값이 중요함\n",
"\n",
"### DIY\n",
"Newton-Raphson method 함수를 구성하시오. 아래 판별 기준을 적용하라.\n",
"- $f(x_i)$ 값이 tolerance 보다 작으면 멈춘다.\n",
"- 근사 상대 오차 $\\epsilon_a$ 가 tolerance 보다 작으면 멈춘다.\n",
"- $f'(x_i) = 0$ 인 경우 에러를 발생시킨다 (`raise ValueError`)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "e63424a2",
"metadata": {},
"outputs": [],
"source": [
"def newton(f, fp, x0, rtol=1e-6):\n",
" # DIY\n",
" ..."
]
},
{
"cell_type": "markdown",
"id": "418a65e8",
"metadata": {},
"source": [
"## Secant Method\n",
"Newton-Raphson 방법은 미분 함수 $f'(x)$ 가 필요하다. 미분을 구하기 힘들 경우 수치적으로 근사한다.\n",
"\n",
":::{figure-md} secant\n",
"
\n",
"\n",
"Secant Method (from Wikipedia)\n",
":::\n",
"\n",
"그림과 같이 ${i-1}$, ${i}$ 번째 결과로 부터 기울기를 대신 사용한다.\n",
"\n",
"$$\n",
"f'(x_i) \\approx \\frac{f(x_i) - f(x_{i-1})}{x_{i} - x_{i-1}}\n",
"$$\n",
"\n",
"이 근사 미분값을 사용하는 방법이 Secant method 이다.\n",
"\n",
"$$\n",
"x_{i+1} = x_{i} - \\frac{f(x_i)}{f'(x_i)} = x_i - \\frac{f(x_i)(x_{i} - x_{i-1})}{f(x_i) - f(x_{i-1})}\n",
"$$\n",
"\n",
"### 예제 1\n",
"다음 함수 $f(x) = e^{-x} - x$ 의 근을 구하시오. Secant Method로 해석해보자.\n",
"(초기 값은 $x_{-1} = 0, x_0 = 1$ 로 한다.)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "3ea29bcc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 0.6127, err = -3.873e-01\n",
"x1 = 0.5638, err = -4.886e-02\n",
"x1 = 0.5672, err = 3.332e-03\n",
"x1 = 0.5671, err = -2.705e-05\n",
"x1 = 0.5671, err = -1.620e-08\n"
]
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : np.exp(-x) - x\n",
"\n",
"xm, x0 = 0.0, 1.0\n",
"x = x0\n",
"\n",
"Err = []\n",
"itmax = 5\n",
"for i in range(itmax):\n",
" # Scecant Method\n",
" xn = x - f(x)*(x - xm) / (f(x) - f(xm))\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" xm = x\n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "markdown",
"id": "f9dd57f6",
"metadata": {},
"source": [
"### 예제 2\n",
"다음 함수 $f(x) = log(x)$ 의 근을 구하시오. Secant Method로 해석해보자.\n",
"(초기 값은 $x_{-1} = 0.5, x_0 = 5$ 로 한다.)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "40f2c61f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1 = 1.8546, err = -3.145e+00\n",
"x1 = -0.1044, err = -1.959e+00\n",
"x1 = nan, err = nan\n",
"x1 = nan, err = nan\n",
"x1 = nan, err = nan\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/tmp/ipykernel_382/3825093768.py:2: RuntimeWarning: invalid value encountered in log\n",
" f = lambda x : np.log(x)\n"
]
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : np.log(x)\n",
"\n",
"xm, x0 = 0.5, 5\n",
"x = x0\n",
"\n",
"Err = []\n",
"itmax = 5\n",
"for i in range(itmax):\n",
" # Scecant Method\n",
" xn = x - f(x)*(x - xm) / (f(x) - f(xm))\n",
" \n",
" # Compute difference\n",
" Eps = xn - x\n",
" \n",
" xm = x\n",
" x = xn\n",
" print(\"x1 = {:.4f}, err = {:.3e}\".format(x, Eps))"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "6c3a71e1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" converged: True\n",
" flag: converged\n",
" function_calls: 44\n",
" iterations: 42\n",
" root: 0.9999999999998863\n",
" method: bisect"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"root_scalar(f, bracket=[0.5, 5], method='bisect')"
]
},
{
"cell_type": "markdown",
"id": "8f611118",
"metadata": {},
"source": [
"### 특징\n",
"- 미분 함수 없이 계산 가능하다.\n",
"- 수렴 속도가 Newton-Raphson과 거의 비슷하나 약간 늦다.\n",
"- 발산할 수 있다."
]
},
{
"cell_type": "markdown",
"id": "eee2ef8f",
"metadata": {},
"source": [
"### DIY\n",
"Scecant method 함수를 구성하시오. Newton-Raphson 과 동일한 판별 조건을 사용하자."
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "ac943e8a",
"metadata": {},
"outputs": [],
"source": [
"def secant(f, x0, x1, rtol=1e-6):\n",
" # DIY\n",
" ..."
]
},
{
"cell_type": "markdown",
"id": "d4e982bb",
"metadata": {},
"source": [
"## SciPy 활용"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "3df0379f",
"metadata": {},
"outputs": [],
"source": [
"# root_scalar?"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "a7dd9b3a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" converged: True\n",
" flag: converged\n",
" function_calls: 8\n",
" iterations: 4\n",
" root: 0.5671432904097838\n",
" method: newton"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function and derivatives\n",
"f = lambda x : np.exp(-x) - x\n",
"fp = lambda x: -np.exp(-x) - 1\n",
"\n",
"root_scalar(f, fprime=fp, x0=1, method='newton')"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "86732f69",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" converged: True\n",
" flag: converged\n",
" function_calls: 7\n",
" iterations: 6\n",
" root: 0.567143290409784\n",
" method: secant"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"root_scalar(f, x0=0, x1=1, method='secant')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}